| 1 | #include "machdefs.h"
 | 
|---|
| 2 | 
 | 
|---|
| 3 | #include <math.h>
 | 
|---|
| 4 | #include <iostream>
 | 
|---|
| 5 | #include <fstream>
 | 
|---|
| 6 | 
 | 
|---|
| 7 | #include "tarrinit.h"
 | 
|---|
| 8 | #include "array.h"
 | 
|---|
| 9 | 
 | 
|---|
| 10 | #include "pexceptions.h"
 | 
|---|
| 11 | #include "srandgen.h"
 | 
|---|
| 12 | 
 | 
|---|
| 13 | #include "sophyainit.h"
 | 
|---|
| 14 | #include "ctimer.h"
 | 
|---|
| 15 | #include "timing.h"
 | 
|---|
| 16 | 
 | 
|---|
| 17 | #include "smtx.h"
 | 
|---|
| 18 | 
 | 
|---|
| 19 | //--------------------------------------------------------
 | 
|---|
| 20 | //   Test d'expression template sur tableaux 
 | 
|---|
| 21 | //       Reza - Avril 2003
 | 
|---|
| 22 | //--------------------------------------------------------
 | 
|---|
| 23 | 
 | 
|---|
| 24 | //-------  Le code de smtx.cc
 | 
|---|
| 25 | 
 | 
|---|
| 26 | #include "smtx.cc"
 | 
|---|
| 27 | 
 | 
|---|
| 28 | //-------
 | 
|---|
| 29 | // ------------------------------------------------
 | 
|---|
| 30 | //   programme de test template expressions 
 | 
|---|
| 31 | // ------------------------------------------------
 | 
|---|
| 32 | 
 | 
|---|
| 33 | int main(int narg, char* arg[])
 | 
|---|
| 34 | {
 | 
|---|
| 35 | 
 | 
|---|
| 36 |   SophyaInit();
 | 
|---|
| 37 |   InitTim();   // Initializing the CPU timer
 | 
|---|
| 38 | 
 | 
|---|
| 39 |   if (narg < 3) {
 | 
|---|
| 40 |     cout << "\n te/ TemplateExpression Test - missing arguments\n" 
 | 
|---|
| 41 |          << "     Usage: te NLoop NRow NCol \n"  << endl;
 | 
|---|
| 42 |     return 1;
 | 
|---|
| 43 |   }
 | 
|---|
| 44 |   int N = 10;
 | 
|---|
| 45 |   int nrow = 2; 
 | 
|---|
| 46 |   int ncol = 2;
 | 
|---|
| 47 | 
 | 
|---|
| 48 |   N = atoi(arg[1]);
 | 
|---|
| 49 | 
 | 
|---|
| 50 |   nrow = atoi(arg[2]);
 | 
|---|
| 51 |   ncol = atoi(arg[3]);
 | 
|---|
| 52 | 
 | 
|---|
| 53 |   cout << " te/ TemplateExpression Test - NLoop = " << N 
 | 
|---|
| 54 |        << "  NRow=" << nrow << " NCol= " << ncol 
 | 
|---|
| 55 |        << "  Size=nr*nc= " << nrow*ncol << endl;
 | 
|---|
| 56 | 
 | 
|---|
| 57 |   try {
 | 
|---|
| 58 |     cout << " ---- Verification operation ------ " << endl;
 | 
|---|
| 59 |     SimpleMatrix<int_4> mx1(4,6);
 | 
|---|
| 60 |     mx1 = 3;
 | 
|---|
| 61 |     mx1 *= 2;
 | 
|---|
| 62 |     cout << " mx1=3; mx1*=2; : " << endl;
 | 
|---|
| 63 |     cout << mx1 << endl;
 | 
|---|
| 64 |     SimpleMatrix<int_4> mx2(4,6);
 | 
|---|
| 65 |     mx2.Init(1,1);
 | 
|---|
| 66 |     SimpleMatrix<int_4> mx3(4,6), mx4;
 | 
|---|
| 67 |     mx3.Init(1000,100);
 | 
|---|
| 68 |     mx4 = mx3+mx2*2 ;
 | 
|---|
| 69 |     cout << " mx2, mx3, mx4=mx3+2mx2 " << endl;
 | 
|---|
| 70 |     cout << mx2 << endl << mx3 << endl << mx4 << endl;
 | 
|---|
| 71 | 
 | 
|---|
| 72 |     cout << ">> Test with SimpleMatrix<T>/JET " << endl;
 | 
|---|
| 73 |     PrtTim(" (0) -- ");
 | 
|---|
| 74 | 
 | 
|---|
| 75 |     int i,j,k;
 | 
|---|
| 76 |     {
 | 
|---|
| 77 |     SimpleMatrix<r_8> m1(nrow, ncol);
 | 
|---|
| 78 |     SimpleMatrix<r_8> m2(nrow, ncol);
 | 
|---|
| 79 |     SimpleMatrix<r_8> m3(nrow, ncol);
 | 
|---|
| 80 |     SimpleMatrix<r_8> m4(nrow, ncol);
 | 
|---|
| 81 |     SimpleMatrix<r_8> m5(nrow, ncol);
 | 
|---|
| 82 |     for(k=0; k<N; k++) {
 | 
|---|
| 83 |       double xxr = frand01();
 | 
|---|
| 84 |       double yyr = 2.*frand01();
 | 
|---|
| 85 |       for(i=0; i<nrow; i++)
 | 
|---|
| 86 |         for(j=0; j<ncol; j++) {
 | 
|---|
| 87 |           m1(i,j) = k*300+10.*i+j+xxr;  
 | 
|---|
| 88 |           m2(i,j) = k*550+20.*i+2.*j+yyr;
 | 
|---|
| 89 |           m3(i,j) = k*860.+40.*i+7.*j+yyr*3.14;
 | 
|---|
| 90 |         }
 | 
|---|
| 91 |     }
 | 
|---|
| 92 | 
 | 
|---|
| 93 |     cout << " (1) Element access done ---" << endl;
 | 
|---|
| 94 |     PrtTim(" (1) -- ");
 | 
|---|
| 95 |     // Calcul m1*c1 + m2*c2 + m3*c3;
 | 
|---|
| 96 |     for(k=0; k<N; k++) {
 | 
|---|
| 97 |       double c1 = frand01() + 1.2;
 | 
|---|
| 98 |       double c2 = frand01() + 3.5;
 | 
|---|
| 99 |       double c3 = frand01() + 6.7;
 | 
|---|
| 100 |       m5 = m1;  m5.MulCst(c1); 
 | 
|---|
| 101 |       m4 = m2;  m4.MulCst(c2);
 | 
|---|
| 102 |       m5.AddElt(m4);
 | 
|---|
| 103 |       m4 = m3;  m3.MulCst(c2);
 | 
|---|
| 104 |       m5.AddElt(m4);
 | 
|---|
| 105 |     }
 | 
|---|
| 106 |     cout << " (2) Add/MulElt/Cst m1*c1 + m2*c2 + m3*c3 done ---" << endl;
 | 
|---|
| 107 |     PrtTim(" (2) -- ");
 | 
|---|
| 108 | 
 | 
|---|
| 109 |     for(k=0; k<N; k++) {
 | 
|---|
| 110 |       double c1 = frand01() + 1.2;
 | 
|---|
| 111 |       double c2 = frand01() + 3.5;
 | 
|---|
| 112 |       double c3 = frand01() + 6.7;
 | 
|---|
| 113 |       m5 = m1*c1+m2*c2+m3*c3;
 | 
|---|
| 114 |     }
 | 
|---|
| 115 |     cout << " (3) JET: m1*c1 + m2*c2 + m3*c3 done ---" << endl;
 | 
|---|
| 116 |     PrtTim(" (3) -- ");
 | 
|---|
| 117 |     // Calcul m1*c1 + m1*(m2+c2) + m2*c3; 
 | 
|---|
| 118 |     for(k=0; k<N; k++) {
 | 
|---|
| 119 |       double c1 = frand01() + 1.2;
 | 
|---|
| 120 |       double c2 = frand01() + 3.5;
 | 
|---|
| 121 |       double c3 = frand01() + 6.7;
 | 
|---|
| 122 |       m3 = m1;  m3.MulCst(c1); 
 | 
|---|
| 123 |       m4 = m2;  m4.AddCst(c2); m4.MulElt(m1);
 | 
|---|
| 124 |       m5 = m2;  m5.MulCst(c3);
 | 
|---|
| 125 |       m5.AddElt(m3); m5.AddElt(m4);
 | 
|---|
| 126 |     }
 | 
|---|
| 127 |     cout << " (4) Add/MulElt/Cst m1*c1 + m1*(m2+c2) + m2*c3 done ---" << endl;
 | 
|---|
| 128 |     PrtTim(" (4) --");
 | 
|---|
| 129 | 
 | 
|---|
| 130 |     for(k=0; k<N; k++) {
 | 
|---|
| 131 |       double c1 = frand01() + 1.2;
 | 
|---|
| 132 |       double c2 = frand01() + 3.5;
 | 
|---|
| 133 |       double c3 = frand01() + 6.7;
 | 
|---|
| 134 |       //      m4 = SMExprAdd< r_8, SMExprMtx<r_8>, SMExprMtx<r_8> >( SMExprMtx<r_8>(m1), SMExprMtx<r_8> (m2) );  
 | 
|---|
| 135 |     // Calcul m1*c1 + m1*(m2+c2) + m2*c3; 
 | 
|---|
| 136 |       m5 = m1*c1 + m1*(m2+c2) + m2*c3;
 | 
|---|
| 137 |     }
 | 
|---|
| 138 |     cout << " (5) JET: m5 = m1*c1 + m1*(m2+c2) + m2*c3 ---" << endl;
 | 
|---|
| 139 |     PrtTim(" (5) -- ");
 | 
|---|
| 140 | 
 | 
|---|
| 141 |     for(k=0; k<N; k++) {
 | 
|---|
| 142 |       double c1 = frand01() + 1.2;
 | 
|---|
| 143 |       double c2 = frand01() + 3.5;
 | 
|---|
| 144 |       double c3 = frand01() + 6.7;
 | 
|---|
| 145 |     // Calcul m1*c1 + m1*(m2+c2) + m2*c3; 
 | 
|---|
| 146 |       m5 = m1*c1 + m2 + c2*Sin(c3*m3+m1);
 | 
|---|
| 147 |     }
 | 
|---|
| 148 |     cout << " (6) JET: m5 = m1*c1 + m2 + c2*Sin(c3*m3+m1) ---" << endl;
 | 
|---|
| 149 |     PrtTim(" (6) -- ");
 | 
|---|
| 150 | 
 | 
|---|
| 151 |     }
 | 
|---|
| 152 | 
 | 
|---|
| 153 |     
 | 
|---|
| 154 |     // Test avec les TArray / TMatrix
 | 
|---|
| 155 |     {
 | 
|---|
| 156 |     cout << " >>>>> test with SOPHYA::TArray / TMatrix " << endl;
 | 
|---|
| 157 |     TMatrix<r_8> m1(nrow, ncol);
 | 
|---|
| 158 |     TMatrix<r_8> m2(nrow, ncol);
 | 
|---|
| 159 |     TMatrix<r_8> m3(nrow, ncol);
 | 
|---|
| 160 |     TMatrix<r_8> m4(nrow, ncol);
 | 
|---|
| 161 |     TMatrix<r_8> m5(nrow, ncol);
 | 
|---|
| 162 |     for(k=0; k<N; k++) {
 | 
|---|
| 163 |       double xxr = frand01();
 | 
|---|
| 164 |       double yyr = 2.*frand01();
 | 
|---|
| 165 |       for(i=0; i<nrow; i++)
 | 
|---|
| 166 |         for(j=0; j<ncol; j++) {
 | 
|---|
| 167 |           m1(i,j) = k*300+10.*i+j+xxr;  
 | 
|---|
| 168 |           m2(i,j) = k*550+20.*i+2.*j+yyr;
 | 
|---|
| 169 |           m3(i,j) = k*860.+40.*i+7.*j+yyr*3.14;
 | 
|---|
| 170 |         }
 | 
|---|
| 171 |     }
 | 
|---|
| 172 |     cout << " (21) Element access done ---" << endl;
 | 
|---|
| 173 |     PrtTim(" (21) -- ");
 | 
|---|
| 174 |     // Calcul m1*c1 + m2*c2 + m3*c3;
 | 
|---|
| 175 |     for(k=0; k<N; k++) {
 | 
|---|
| 176 |       double c1 = frand01() + 1.2;
 | 
|---|
| 177 |       double c2 = frand01() + 3.5;
 | 
|---|
| 178 |       double c3 = frand01() + 6.7;
 | 
|---|
| 179 |       m5 = m1;  m5.Mul(c1); 
 | 
|---|
| 180 |       m4 = m2;  m4.Mul(c2);
 | 
|---|
| 181 |       m5.AddElt(m4);
 | 
|---|
| 182 |       m4 = m3;  m3.Mul(c2);
 | 
|---|
| 183 |       m5.AddElt(m4);
 | 
|---|
| 184 |     }
 | 
|---|
| 185 |     cout << " (22) Add/MulElt/Cst m1*c1 + m2*c2 + m3*c3 done ---" << endl;
 | 
|---|
| 186 |     PrtTim(" (22) -- ");
 | 
|---|
| 187 | 
 | 
|---|
| 188 |     for(k=0; k<N; k++) {
 | 
|---|
| 189 |       double c1 = frand01() + 1.2;
 | 
|---|
| 190 |       double c2 = frand01() + 3.5;
 | 
|---|
| 191 |       double c3 = frand01() + 6.7;
 | 
|---|
| 192 |       m5 = m1*c1+m2*c2+m3*c3;
 | 
|---|
| 193 |     }
 | 
|---|
| 194 |     cout << " (23) m1*c1 + m2*c2 + m3*c3 done ---" << endl;
 | 
|---|
| 195 |     PrtTim(" (23) -- ");
 | 
|---|
| 196 | 
 | 
|---|
| 197 |     for(k=0; k<N; k++) {
 | 
|---|
| 198 |       double c1 = frand01() + 1.2;
 | 
|---|
| 199 |       double c2 = frand01() + 3.5;
 | 
|---|
| 200 |       double c3 = frand01() + 6.7;
 | 
|---|
| 201 |       m5 = m1*c1 + m2 + c2*Sin(c3*m3+m1);
 | 
|---|
| 202 |     }
 | 
|---|
| 203 |     cout << " (26) m5 = m1*c1 + m2 + c2*Sin(c3*m3+m1) ---" << endl;
 | 
|---|
| 204 |     PrtTim(" (26) -- ");
 | 
|---|
| 205 | 
 | 
|---|
| 206 |     }
 | 
|---|
| 207 | 
 | 
|---|
| 208 |     // Test avec des double * x;
 | 
|---|
| 209 |     {
 | 
|---|
| 210 |     cout << " >>>>> test with r_8 * p = new double[size] " << endl;
 | 
|---|
| 211 |     r_8 * m1 = new r_8[nrow*ncol];
 | 
|---|
| 212 |     r_8 * m2 = new r_8[nrow*ncol];
 | 
|---|
| 213 |     r_8 * m3 = new r_8[nrow*ncol];
 | 
|---|
| 214 |     r_8 * m5 = new r_8[nrow*ncol];
 | 
|---|
| 215 |     for(k=0; k<N; k++) {
 | 
|---|
| 216 |       double xxr = frand01();
 | 
|---|
| 217 |       double yyr = 2.*frand01();
 | 
|---|
| 218 |       for(i=0; i<nrow*ncol; i++) {
 | 
|---|
| 219 |           m1[i] = k*300+10.*i+j+xxr;    
 | 
|---|
| 220 |           m2[i] = k*550+20.*i+2.*j+yyr;
 | 
|---|
| 221 |           m3[i] = k*860.+40.*i+7.*j+yyr*3.14;
 | 
|---|
| 222 |       }
 | 
|---|
| 223 |     }
 | 
|---|
| 224 |     cout << " (31) p[i] Element access done ---" << endl;
 | 
|---|
| 225 |     PrtTim(" (31) -- ");
 | 
|---|
| 226 |     // Calcul m1*c1 + m2*c2 + m3*c3;
 | 
|---|
| 227 |     for(k=0; k<N; k++) {
 | 
|---|
| 228 |       double c1 = frand01() + 1.2;
 | 
|---|
| 229 |       double c2 = frand01() + 3.5;
 | 
|---|
| 230 |       double c3 = frand01() + 6.7;
 | 
|---|
| 231 |       for(i=0; i<nrow*ncol; i++) 
 | 
|---|
| 232 |         m5[i] = c1*m1[i]+c2*m2[i]+c3*m3[i];
 | 
|---|
| 233 |     }
 | 
|---|
| 234 |     cout << " (32) p[i] m1*c1 + m2*c2 + m3*c3 done ---" << endl;
 | 
|---|
| 235 |     PrtTim(" (32) -- ");
 | 
|---|
| 236 | 
 | 
|---|
| 237 |     delete[] m1;
 | 
|---|
| 238 |     delete[] m2;
 | 
|---|
| 239 |     delete[] m3;
 | 
|---|
| 240 |     delete[] m5;
 | 
|---|
| 241 |     }
 | 
|---|
| 242 |   }
 | 
|---|
| 243 |   catch (PThrowable exc) {
 | 
|---|
| 244 |     cerr << " te/catched Exception " << exc.Msg() << endl;
 | 
|---|
| 245 |   }  
 | 
|---|
| 246 |   catch (...) {
 | 
|---|
| 247 |     cerr << " te/catched unknown (...) exception " << endl; 
 | 
|---|
| 248 |   }  
 | 
|---|
| 249 | 
 | 
|---|
| 250 |   cout << "\n --------------------------------------------------------" << endl;
 | 
|---|
| 251 |   PrtTim("--- End of te ---");
 | 
|---|
| 252 |   cout << " ----  END of te/TemplateExpression programme ------ " << endl;
 | 
|---|
| 253 |   return 0;
 | 
|---|
| 254 | }
 | 
|---|
| 255 |   
 | 
|---|