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 }