1 double f(double& d)
 2 {
 3     d*=2;
 4     return d;
 5 }
 6
 7 void f1(double&d)
 8 {
 9     d*=2;
10 }
11
12 double f2(double& d)
13 {
14     d*=2;
15     return d;
16 }

1 class f3{
2     int i;
3 public:
4     f3(int i1):i(i1){}
5     //double operator()(double d){return d*i;}//三者不能共存
6     //void operator()(double&d){d*=i;}
7     double operator()(double&d){d*=i;return d;}
8 };

  1 #include <iostream>
  2 #include <vector>
  3 #include <string>
  4 #include <list>
  5 #include<fstream>
  6 #include <set>
  7 #include<algorithm>
  8 #include<stdexcept>
  9 #include <map>
 10 #include<Matrix.h>
 11 #include<MatrixIO.h>
 12 using namespace Numeric_lib;
 13 using namespace std;
 14
 15 typedef Matrix<double,2>Matrix1;
 16 typedef Matrix<double,1>vector1;
 17
 18 void elim_with_partial_pivot(Matrix1&A,vector1&b)
 19 {
 20     const Index n=A.dim1();
 21
 22     for(Index j=0;j<n;++j)
 23     {
 24         Index pivot_row=j;
 25         for(Index k=j+1;k<n;++k)
 26             if(abs(A(k,j))>abs(A(pivot_row,j)))
 27                 pivot_row=j;
 28
 29         if(pivot_row!=j)
 30         {
 31             A.swap_rows(j,pivot_row);
 32             std::swap(b(j),b(pivot_row));
 33         }
 34         for(Index i=j+1;i<n;++i)
 35         {
 36             const double pivot=A(j,j);
 37             if(pivot==0)
 38                 cerr<<"hehe"<<endl;
 39             const double mult=A(i,j)/pivot;
 40             A[i].slice(j)=scale_and_add(A[j].slice(j),-mult,A[i].slice(j));
 41             b(i)-=mult*b(j);
 42         }
 43
 44     }
 45 }
 46
 47 void classical_elimination(Matrix1&A,vector1&b)
 48 {
 49     const Index n=A.dim1();
 50     for(Index j=0;j<n-1;++j)
 51     {
 52         const double pivot=A(j,j);
 53         if(pivot==0)
 54             throw exception();
 55
 56         for(Index i=j+1;i<n;++i)
 57         {
 58             const double mult =A(i,j)/pivot;
 59             A[i].slice(j)=scale_and_add(A[j].slice(j),-mult,A[i].slice(j));
 60             b(i)-=mult*b(j);
 61         }
 62     }
 63 }
 64
 65 vector1 back_substitution(const Matrix1&A,const vector1&b)
 66 {
 67     const Index n=A.dim1();
 68     vector1 x(n);
 69     for(Index i=n-1;i>=0;--i)
 70     {
 71         double s=b(i)-dot_product(A[i].slice(i+1),x.slice(i+1));
 72         if(double m=A(i,i))
 73             x(i)=s/m;
 74         else
 75             throw exception();
 76     }
 77     return x;
 78 }
 79
 80 vector1 classical_gaussian_elimination(Matrix1 A,vector1 b)
 81 {
 82     //classical_elimination(A,b);
 83     elim_with_partial_pivot(A,b);
 84     return back_substitution(A,b);
 85 }
 86
 87 vector1 operator*(const Matrix1&m,const vector1&u)
 88 {
 89     const Index n=m.dim1();
 90     vector1 v(n);
 91     for(Index i=0;i<n;++i)
 92         v(i)=dot_product(m[i],u);
 93     return v;
 94 }
 95
 96 void solve_random_system()
 97 {
 98     Matrix1 A(2,2);
 99     A(0,0)=1;
100     A(0,1)=0;
101     A(1,0)=0;
102     A(1,1)=1;
103     vector1 b(2);
104     b(0)=5;b(1)=6;
105
106     cout<<A<<endl<<b<<endl;
107     try{
108         vector1 x=classical_gaussian_elimination(A,b);
109         cout<<endl<<x<<endl;
110         vector1 v=A*x;
111         cout<<endl<<v<<endl;
112     }
113     catch(const exception&e)
114     {
115         cerr<<e.what()<<endl;
116     }
117 }
118
119 int main()
120 {
121     solve_random_system();
122
123     while(1);
124     return 0;
125
126 }

1 template<class T,class P>
2 Matrix<T> apply1(P p,Matrix<T,1> a)
3 {
4     Matrix<T>t1=a;
5     for(Index i=0;i<t1.size();i++)
6         t1(i)=p(t1(i));
7     return t1;
8 }

 1 template<class T>
 2 void swap_columns(Matrix<T,2>&m,Index i, Index j)
 3 {
 4     if (i == j) return;
 5
 6     T t;
 7     int k=m.dim1();
 8     for(int h=0;h<k;h++)
 9     {
10         t=m[h][i];
11         m[h][i]=m[h][j];
12         m[h][j]=t;
13     }
14 }