00001 #include <iostream.h> 00002 #include<gslwrap/matrix_float.h> 00003 #include<gslwrap/matrix_double.h> 00004 #include<gslwrap/vector_float.h> 00005 #include<gslwrap/vector_int.h> 00006 #include<gslwrap/vector_double.h> 00007 #include<gslwrap/random_generator.h> 00008 using namespace gsl; 00009 00010 void 00011 test_vector() 00012 { 00013 gsl::vector vd(10); 00014 gsl::vector_float vf(10); 00015 gsl::vector_int vi(10); 00016 } 00017 00018 // Testing creation, dimensioning and assigning functions: 00019 void 00020 vector_float_test() 00021 { 00022 vector_float test; 00023 vector_float mal(10); 00024 for (int i=0; i<10; i++) 00025 {mal[i] = i;} 00026 00027 test = mal; 00028 cout << "Test vector should be:" << endl <<mal<< endl; 00029 cout << "Test vector is:" << endl <<test<< endl; 00030 00031 if (mal != test) 00032 cout << "vector_float ERROR Test did not pass !! --------------" << endl; 00033 else 00034 cout << "vector_float Test passed !! --------------" << endl; 00035 } 00036 00037 void 00038 VectorViewTest() 00039 { 00040 cout << "VectorViewTest --------------------------" << endl; 00041 int i,j; 00042 int size=5; 00043 matrix_float m(size, size); 00044 matrix_float ver(size, size);ver.set_all(5.0); 00045 vector_float v(size); 00046 for (i=0;i<size;i++) 00047 { 00048 v[i]=size-(i+1); 00049 for (j=0;j<size;j++) 00050 m(i,j)=(i+1); 00051 } 00052 cout << "m=" << endl << m << endl; 00053 cout << "v=" << endl << v << endl; 00054 00055 vector_float_view col_viewer = m.column(3); 00056 vector_float_view viewer = m.row(3); 00057 cout << "viewer (3.row before adding)=" << endl << viewer << endl; 00058 cout << "col_viewer (3.col before adding)=" << endl << col_viewer << endl; 00059 for (i=0;i<m.get_cols();i++) 00060 m.column(i) += v; 00061 00062 cout << "viewer (3.row after adding) =" << endl << viewer << endl; 00063 cout << "col_viewer (3.col after adding)=" << endl << col_viewer << endl; 00064 cout << "\"m+v\"" << endl << m << endl; 00065 00066 col_viewer.change_view(viewer); 00067 cout << "col_viewer after changing to viewer=" << endl << col_viewer << endl; 00068 00069 if (m!=ver || m.column(0)!=viewer || col_viewer!=viewer) 00070 cout << "VectorViewTest failed !! ----------------" << endl; 00071 else 00072 cout << "VectorViewTest passed !! ----------------" << endl; 00073 00074 } 00075 00076 void VectorViewTest2() 00077 { 00078 cout << "VectorViewTest2 ----------------------------" << endl; 00079 matrix m1(10,10); 00080 matrix m2(10,10); 00081 vector v(10); 00082 int i; 00083 00084 for (i=0;i<10; i++){v[i]=i*2;} 00085 00086 for (i=0;i<10; i++) 00087 { 00088 m1.column(i) = v; 00089 m2.column(i) = v; 00090 } 00091 00092 for (i=0;i<10; i++) 00093 m1.column(i) -= m2.column(i); 00094 00095 cout << "Sum m1-m2 (m1==m2) : " << m1.sum() << endl; 00096 if (m1.sum()) 00097 { 00098 cout << "VectorViewTest2 Failed !! (operator -= ?) " << endl; 00099 exit(-1); 00100 } 00101 else 00102 cout << "VectorViewTest2 Passed !! -------------------------" << endl; 00103 } 00104 00105 void VectorViewTest3() 00106 { 00107 vector v(10); 00108 vector_view view1 = v; 00109 vector_view view2 = v; 00110 view2 = view1; 00111 } 00112 //some test calls to gsl: 00113 void 00114 gsl_function_call_test() 00115 { 00116 int size=5; 00117 matrix m(size,size); 00118 for (int i=0;i<size;i++) 00119 for (int j=0;j<size;j++) 00120 m(i,j) = (i+1)*(j+1); 00121 double lndet = gsl_linalg_LU_lndet(m.gslobj()); 00122 cout << "m = " << m << endl; 00123 cout << "lndet(m) = " << lndet << endl; 00124 00125 gsl::vector sol(size); 00126 for (int k=0;k<size;k++) 00127 sol[k]=k+1; 00128 cout << "Solution of " << m << " x= " << sol <<endl; 00129 00130 gsl::vector tau(size); 00131 gsl_linalg_QR_decomp(m.gslobj(), tau.gslobj()); 00132 gsl_linalg_QR_svx(m.gslobj(), tau.gslobj(), sol.gslobj()); 00133 00134 cout << "sol= " << sol << endl; 00135 } 00136 00137 void 00138 RandomNumberGeneratorTest() 00139 { 00140 cout << "RandomNumberGeneratorTest -------------------" << endl; 00141 random_generator generator; // default constructor takes seed from ENV variable 00142 int nSamples=10000; 00143 gsl::vector samples(nSamples); 00144 for (int i=0;i<nSamples;i++) 00145 samples[i] = generator.uniform(); 00146 00147 double mean=samples.sum()/(float)nSamples; 00148 cout << "Mean of sampling : " << mean << endl; 00149 // should be close to 0.5 00150 if (mean <0.51 && mean >0.49) 00151 cout << "RandomNumberGeneratorTest Passed !! -------------" << endl; 00152 else 00153 { 00154 cout << "RandomNumberGeneratorTest Failed !! -------------" << endl; 00155 exit(-1); 00156 } 00157 } 00158 00159 int main( int argc, char **argv ) 00160 { 00161 matrix a( 3, 3 ); 00162 matrix b( 3, 3 ); 00163 srand48(10); 00164 00165 int i, j; 00166 for(i=0;i<a.size1();i++) 00167 for(j=0;j<a.size2();j++) a(i,j)=drand48(); 00168 for(i=0;i<b.size1();i++) 00169 for(j=0;j<b.size2();j++) b(i,j)=drand48(); 00170 00171 cout << "a:" << endl << a << endl; 00172 cout << "b:" << endl << b << endl; 00173 00174 cout << "a + b:" << endl << a + b << endl; 00175 cout << "a - b:" << endl << a - b << endl; 00176 cout << "a + 1.0:" << endl << a + 1.0 << endl; 00177 cout << "1.0 - a:" << endl << 1.0 - a << endl; 00178 cout << "a * 10.0:" << endl << a * 10.0 << endl; 00179 cout << "10.0 * a:" << endl << 10.0 * a << endl; 00180 cout << "a / 10.0:" << endl << a / 10.0 << endl; 00181 cout << "a * b:" << endl << a * b << endl; 00182 cout << "a.LU_decomp():" << endl << a.LU_decomp() << endl; 00183 cout << "a.LU_invert():" << endl << a.LU_invert() << endl; 00184 cout << "a.LU_invert() * a:" << endl << a.LU_invert() * a << endl; 00185 00186 // call class test functions (there should be one per class) 00187 vector_float_test(); 00188 gsl_function_call_test(); 00189 VectorViewTest(); 00190 RandomNumberGeneratorTest(); 00191 VectorViewTest2(); 00192 VectorViewTest3(); 00193 }