Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cxx_tmpl_main_comp.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Teuchos: Common Tools Package
4 //
5 // Copyright 2004 NTESS and the Teuchos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 // Kris
11 // 07.24.03 -- Initial checkin
12 // 08.08.03 -- All test suites except for TRSM are finished.
13 // 08.14.03 -- The test suite for TRSM is finished (Heidi).
14 /*
15 
16 This test program is intended to check an experimental default type (e.g. mp_real) against an "officialy supported" control type (e.g. double). For each test, the program will generate the appropriate scalars and randomly-sized vectors and matrices of random numbers for the routine being tested. All of the input data for the experimental type is casted into the control type, so in theory both BLAS routines should get the same input data. Upon return, the experimental output data is casted back into the control type, and the results are compared; if they are equal (within a user-definable tolerance) the test is considered successful.
17 
18 The test routine for TRSM is still being developed; all of the others are more or less finalized.
19 
20 */
21 
22 #include "Teuchos_BLAS.hpp"
23 #include "Teuchos_Time.hpp"
24 #include "Teuchos_Version.hpp"
25 #include "Teuchos_GlobalMPISession.hpp"
26 
27 using Teuchos::BLAS;
28 using Teuchos::ScalarTraits;
29 
30 // OType1 and OType2 define the ordinal datatypes for which BLAS output will be compared.
31 // The difference in OType should enable the comparison of the templated routines with the "officially" supported BLAS.
32 
33 // Define the scalar type
34 #ifdef HAVE_TEUCHOS_COMPLEX
35 #define SType std::complex<double>
36 #else
37 #define SType double
38 #endif
39 
40 // Define the ordinal type
41 #define OType1 long int
42 #define OType2 int
43 
44 // MVMIN/MAX define the minimum and maximum dimensions of generated matrices and vectors, respectively.
45 // These are well within the range of OType1 and OType2
46 #define MVMIN 2
47 #define MVMAX 20
48 // SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and std::vector elements and scalars:
49 // random numbers in [-SCALARMAX, SCALARMAX] will be generated.
50 // Set SCALARMAX to a floating-point value (e.g. 10.0) to enable floating-point random number generation, such that
51 // random numbers in (-SCALARMAX - 1, SCALARMAX + 1) will be generated.
52 #ifdef HAVE_TEUCHOS_COMPLEX
53 #define SCALARMAX SType(10,0)
54 #else
55 #define SCALARMAX SType(10)
56 #endif
57 // These define the number of tests to be run for each individual BLAS routine.
58 #define ROTGTESTS 5
59 #define ROTTESTS 5
60 #define ASUMTESTS 5
61 #define AXPYTESTS 5
62 #define COPYTESTS 5
63 #define DOTTESTS 5
64 #define IAMAXTESTS 5
65 #define NRM2TESTS 5
66 #define SCALTESTS 5
67 #define GEMVTESTS 5
68 #define GERTESTS 5
69 #define TRMVTESTS 5
70 #define GEMMTESTS 5
71 #define SYMMTESTS 5
72 #define SYRKTESTS 5
73 #define TRMMTESTS 5
74 #define TRSMTESTS 5
75 
76 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
77 template<typename TYPE>
78 TYPE GetRandom(TYPE, TYPE);
79 
80 // Returns a random integer between the two input parameters, inclusive
81 template<>
82 int GetRandom(int, int);
83 
84 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
85 template<>
86 double GetRandom(double, double);
87 
88 template<typename T>
89 std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
90 
91 template<typename TYPE, typename OTYPE>
92 void PrintVector(TYPE* Vector, OTYPE Size, std::string Name, bool Matlab = 0);
93 
94 template<typename TYPE, typename OTYPE>
95 void PrintMatrix(TYPE* Matrix, OTYPE Rows, OTYPE Columns, OTYPE LDM, std::string Name, bool Matlab = 0);
96 
97 template<typename TYPE>
98 bool CompareScalars(TYPE Scalar1, TYPE Scalar2, typename ScalarTraits<TYPE>::magnitudeType Tolerance );
99 
100 template<typename TYPE, typename OTYPE1, typename OTYPE2>
101 bool CompareVectors(TYPE* Vector1, OTYPE1 Size1, TYPE* Vector2, OTYPE2 Size2, typename ScalarTraits<TYPE>::magnitudeType Tolerance );
102 
103 template<typename TYPE, typename OTYPE1, typename OTYPE2>
104 bool CompareMatrices(TYPE* Matrix1, OTYPE1 Rows1, OTYPE1 Columns1, OTYPE1 LDM1,
105  TYPE* Matrix2, OTYPE2 Rows2, OTYPE2 Columns2, OTYPE2 LDM2,
106  typename ScalarTraits<TYPE>::magnitudeType Tolerance );
107 
108 // Use this to convert the larger ordinal type to the smaller one (nothing inherently makes sure of this).
109 template<typename OTYPE1, typename OTYPE2>
110 OTYPE2 ConvertType(OTYPE1 T1, OTYPE2 T2)
111 {
112  return static_cast<OTYPE2>(T1);
113 }
114 
115 // These functions return a random character appropriate for the BLAS arguments that share their names (uses GetRandom())
116 Teuchos::ESide RandomSIDE();
117 Teuchos::EUplo RandomUPLO();
118 Teuchos::ETransp RandomTRANS();
119 Teuchos::EDiag RandomDIAG();
120 
121 int main(int argc, char *argv[])
122 {
123  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
124  bool verbose = 0;
125  bool debug = 0;
126  bool matlab = 0;
127  bool InvalidCmdLineArgs = 0;
128  int i;
129  OType1 j1, k1;
130  OType2 j2, k2;
131  for(i = 1; i < argc; i++)
132  {
133  if(argv[i][0] == '-')
134  {
135  switch(argv[i][1])
136  {
137  case 'v':
138  if(!verbose)
139  {
140  verbose = 1;
141  }
142  else
143  {
144  InvalidCmdLineArgs = 1;
145  }
146  break;
147  case 'd':
148  if(!debug)
149  {
150  debug = 1;
151  }
152  else
153  {
154  InvalidCmdLineArgs = 1;
155  }
156  break;
157  case 'm':
158  if(!matlab)
159  {
160  matlab = 1;
161  }
162  else
163  {
164  InvalidCmdLineArgs = 1;
165  }
166  break;
167  default:
168  InvalidCmdLineArgs = 1;
169  break;
170  }
171  }
172  }
173 
174  if (verbose)
175  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
176 
177  if(InvalidCmdLineArgs || (argc > 4))
178  {
179  std::cout << "Invalid command line arguments detected. Use the following flags:" << std::endl
180  << "\t -v enables verbose mode (reports number of failed/successful tests)" << std::endl
181  << "\t -d enables debug mode (same as verbose with output of each test, not recommended for large numbers of tests)" << std::endl
182  << "\t -m enables matlab-style output; only has an effect if debug mode is enabled" << std::endl;
183  return 1;
184  }
185  typedef ScalarTraits<SType>::magnitudeType MType;
186  BLAS<OType1, SType> OType1BLAS;
187  BLAS<OType2, SType> OType2BLAS;
188  SType STypezero = ScalarTraits<SType>::zero();
189  SType STypeone = ScalarTraits<SType>::one();
190  SType OType1alpha, OType1beta;
191  SType OType2alpha, OType2beta;
192  SType *OType1A, *OType1B, *OType1C, *OType1x, *OType1y;
193  SType *OType2A, *OType2B, *OType2C, *OType2x, *OType2y;
194  SType OType1ASUMresult, OType1DOTresult, OType1NRM2result, OType1SINresult;
195  SType OType2ASUMresult, OType2DOTresult, OType2NRM2result, OType2SINresult;
196  MType OType1COSresult, OType2COSresult;
197  OType1 incx1, incy1;
198  OType2 incx2, incy2;
199  OType1 OType1IAMAXresult;
200  OType2 OType2IAMAXresult;
201  OType1 TotalTestCount = 1, GoodTestSubcount, GoodTestCount = 0, M1, N1, P1, K1, LDA1, LDB1, LDC1, Mx1, My1;
202  OType2 M2, N2, P2, K2, LDA2, LDB2, LDC2, Mx2, My2;
203  Teuchos::EUplo UPLO;
204  Teuchos::ESide SIDE;
205  Teuchos::ETransp TRANS, TRANSA, TRANSB;
206  Teuchos::EDiag DIAG;
207  MType TOL = 1e-8*ScalarTraits<MType>::one();
208 
209  std::srand(time(NULL));
210 
211  //--------------------------------------------------------------------------------
212  // BEGIN LEVEL 1 BLAS TESTS
213  //--------------------------------------------------------------------------------
214  // Begin ROTG Tests
215  //--------------------------------------------------------------------------------
216  GoodTestSubcount = 0;
217  for(i = 0; i < ROTGTESTS; i++)
218  {
219  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
220  OType2alpha = OType1alpha;
221  OType1beta = GetRandom(-SCALARMAX, SCALARMAX);
222  OType2beta = OType1beta;
223  OType1COSresult = ScalarTraits<MType>::zero();
224  OType2COSresult = OType1COSresult;
225  OType1SINresult = ScalarTraits<SType>::zero();
226  OType2SINresult = OType1SINresult;
227 
228  if(debug)
229  {
230  std::cout << "Test #" << TotalTestCount << " -- ROTG -- " << std::endl;
231  std::cout << "OType1alpha = " << OType1alpha << std::endl;
232  std::cout << "OType2alpha = " << OType2alpha << std::endl;
233  std::cout << "OType1beta = " << OType1beta << std::endl;
234  std::cout << "OType2beta = " << OType2beta << std::endl;
235  }
236  TotalTestCount++;
237  OType1BLAS.ROTG(&OType1alpha, &OType1beta, &OType1COSresult, &OType1SINresult);
238  OType2BLAS.ROTG(&OType2alpha, &OType2beta, &OType2COSresult, &OType2SINresult);
239  if(debug)
240  {
241  std::cout << "OType1 ROTG COS result: " << OType1COSresult << std::endl;
242  std::cout << "OType2 ROTG COS result: " << OType2COSresult << std::endl;
243  std::cout << "OType1 ROTG SIN result: " << OType1SINresult << std::endl;
244  std::cout << "OType2 ROTG SIN result: " << OType2SINresult << std::endl;
245  }
246  if ( !CompareScalars(OType1COSresult, OType2COSresult, TOL) || !CompareScalars(OType1SINresult, OType2SINresult, TOL) )
247  std::cout << "FAILED TEST!!!!!!" << std::endl;
248  GoodTestSubcount += ( CompareScalars(OType1COSresult, OType2COSresult, TOL) &&
249  CompareScalars(OType1SINresult, OType2SINresult, TOL) );
250  }
251  GoodTestCount += GoodTestSubcount;
252  if(verbose || debug) std::cout << "ROTG: " << GoodTestSubcount << " of " << ROTGTESTS << " tests were successful." << std::endl;
253  if(debug) std::cout << std::endl;
254  //--------------------------------------------------------------------------------
255  // End ROTG Tests
256  //--------------------------------------------------------------------------------
257 
258  //--------------------------------------------------------------------------------
259  // Begin ROT Tests
260  //--------------------------------------------------------------------------------
261  GoodTestSubcount = 0;
262  for(i = 0; i < ROTTESTS; i++)
263  {
264  incx1 = GetRandom(-5,5);
265  incy1 = GetRandom(-5,5);
266  if (incx1 == 0) incx1 = 1;
267  if (incy1 == 0) incy1 = 1;
268  incx2 = ConvertType( incx1, incx2 );
269  incy2 = ConvertType( incy1, incy2 );
270  M1 = GetRandom(MVMIN, MVMIN+8);
271  M2 = ConvertType( M1, M2 );
272  Mx1 = M1*std::abs(incx1);
273  My1 = M1*std::abs(incy1);
274  if (Mx1 == 0) { Mx1 = 1; }
275  if (My1 == 0) { My1 = 1; }
276  Mx2 = ConvertType( Mx1, Mx2 );
277  My2 = ConvertType( My1, My2 );
278  OType1x = new SType[Mx1];
279  OType1y = new SType[My1];
280  OType2x = new SType[Mx2];
281  OType2y = new SType[My2];
282  for(j1 = 0, j2 = 0; j1 < Mx1; j1++, j2++)
283  {
284  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
285  OType2x[j2] = OType1x[j1];
286  }
287  for(j1 = 0, j2 = 0; j1 < My1; j1++, j2++)
288  {
289  OType1y[j1] = GetRandom(-SCALARMAX, SCALARMAX);
290  OType2y[j2] = OType1y[j1];
291  }
292  MType c1 = cos(ScalarTraits<SType>::magnitude(GetRandom(-SCALARMAX,SCALARMAX)));
293  MType c2 = c1;
294  SType s1 = sin(ScalarTraits<SType>::magnitude(GetRandom(-SCALARMAX,SCALARMAX)));
295  SType s2 = s1;
296  if(debug)
297  {
298  std::cout << "Test #" << TotalTestCount << " -- ROT -- " << std::endl;
299  std::cout << "c1 = " << c1 << ", s1 = " << s1 << std::endl;
300  std::cout << "c2 = " << c2 << ", s2 = " << s2 << std::endl;
301  std::cout << "incx1 = " << incx1 << ", incy1 = " << incy1 << std::endl;
302  std::cout << "incx2 = " << incx2 << ", incy2 = " << incy2 << std::endl;
303  PrintVector(OType1x, Mx1, "OType1x", matlab);
304  PrintVector(OType1y, My1, "OType1y_before_operation", matlab);
305  PrintVector(OType2x, Mx2, "OType2x", matlab);
306  PrintVector(OType2y, My2, "OType2y_before_operation", matlab);
307  }
308  TotalTestCount++;
309  OType1BLAS.ROT(M1, OType1x, incx1, OType1y, incy1, &c1, &s1);
310  OType2BLAS.ROT(M2, OType2x, incx2, OType2y, incy2, &c2, &s2);
311  if(debug)
312  {
313  PrintVector(OType1y, My1, "OType1y_after_operation", matlab);
314  PrintVector(OType2y, My2, "OType2y_after_operation", matlab);
315  }
316  if ( !CompareVectors(OType1x, Mx1, OType2x, Mx2, TOL) || !CompareVectors(OType1y, My1, OType2y, My2, TOL) )
317  std::cout << "FAILED TEST!!!!!!" << std::endl;
318  GoodTestSubcount += ( CompareVectors(OType1x, Mx1, OType2x, Mx2, TOL) &&
319  CompareVectors(OType1y, My1, OType2y, My2, TOL) );
320  delete [] OType1x;
321  delete [] OType1y;
322  delete [] OType2x;
323  delete [] OType2y;
324  }
325  GoodTestCount += GoodTestSubcount;
326  if(verbose || debug) std::cout << "ROT: " << GoodTestSubcount << " of " << ROTTESTS << " tests were successful." << std::endl;
327  if(debug) std::cout << std::endl;
328  //--------------------------------------------------------------------------------
329  // End ROT Tests
330  //--------------------------------------------------------------------------------
331 
332  //--------------------------------------------------------------------------------
333  // Begin ASUM Tests
334  //--------------------------------------------------------------------------------
335  GoodTestSubcount = 0;
336  ScalarTraits<int>::seedrandom(0);
337  for(i = 0; i < ASUMTESTS; i++)
338  {
339  incx1 = GetRandom(1, MVMAX);
340  incx2 = ConvertType( incx1, incx2 );
341  M1 = GetRandom(MVMIN, MVMAX);
342  M2 = ConvertType( M1, M2 );
343  OType1x = new SType[M1*incx1];
344  OType2x = new SType[M2*incx2];
345  for(j1 = 0, j2 = 0; j2 < M2*incx2; j1++, j2++)
346  {
347  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
348  OType2x[j2] = OType1x[j1];
349  }
350  if(debug)
351  {
352  std::cout << "Test #" << TotalTestCount << " -- ASUM -- " << std::endl;
353  std::cout << "incx1 = " << incx1 << "\t" << "incx2 = " << incx2
354  << "\t" << "M1 = " << M1 << "\t" << "M2 = " << M2 << std::endl;
355  PrintVector(OType1x, M1*incx1, "OType1x", matlab);
356  PrintVector(OType2x, M2*incx2, "OType2x", matlab);
357  }
358  TotalTestCount++;
359  OType1ASUMresult = OType1BLAS.ASUM(M1, OType1x, incx1);
360  OType2ASUMresult = OType2BLAS.ASUM(M2, OType2x, incx2);
361  if(debug)
362  {
363  std::cout << "OType1 ASUM result: " << OType1ASUMresult << std::endl;
364  std::cout << "OType2 ASUM result: " << OType2ASUMresult << std::endl;
365  }
366  if (CompareScalars(OType1ASUMresult, OType2ASUMresult, TOL)==0)
367  std::cout << "FAILED TEST!!!!!!" << std::endl;
368  GoodTestSubcount += CompareScalars(OType1ASUMresult, OType2ASUMresult, TOL);
369 
370  delete [] OType1x;
371  delete [] OType2x;
372  }
373  GoodTestCount += GoodTestSubcount;
374  if(verbose || debug) std::cout << "ASUM: " << GoodTestSubcount << " of " << ASUMTESTS << " tests were successful." << std::endl;
375  if(debug) std::cout << std::endl;
376 
377  //--------------------------------------------------------------------------------
378  // End ASUM Tests
379  //--------------------------------------------------------------------------------
380 
381  //--------------------------------------------------------------------------------
382  // Begin AXPY Tests
383  //--------------------------------------------------------------------------------
384  GoodTestSubcount = 0;
385  for(i = 0; i < AXPYTESTS; i++)
386  {
387  incx1 = GetRandom(1, MVMAX);
388  incy1 = GetRandom(1, MVMAX);
389  incx2 = ConvertType( incx1, incx2 );
390  incy2 = ConvertType( incy1, incy2 );
391  M1 = GetRandom(MVMIN, MVMAX);
392  M2 = ConvertType( M1, M2 );
393  Mx1 = M1*std::abs(incx1);
394  My1 = M1*std::abs(incy1);
395  if (Mx1 == 0) { Mx1 = 1; }
396  if (My1 == 0) { My1 = 1; }
397  Mx2 = ConvertType( Mx1, Mx2 );
398  My2 = ConvertType( My1, My2 );
399  OType1x = new SType[Mx1];
400  OType1y = new SType[My1];
401  OType2x = new SType[Mx2];
402  OType2y = new SType[My2];
403  for(j1 = 0, j2 = 0; j1 < Mx1; j1++, j2++)
404  {
405  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
406  OType2x[j2] = OType1x[j1];
407  }
408  for(j1 = 0, j2 = 0; j1 < My1; j1++, j2++)
409  {
410  OType1y[j1] = GetRandom(-SCALARMAX, SCALARMAX);
411  OType2y[j2] = OType1y[j1];
412  }
413  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
414  OType2alpha = OType1alpha;
415  if(debug)
416  {
417  std::cout << "Test #" << TotalTestCount << " -- AXPY -- " << std::endl;
418  std::cout << "OType1alpha = " << OType1alpha << std::endl;
419  std::cout << "OType2alpha = " << OType2alpha << std::endl;
420  PrintVector(OType1x, Mx1, "OType1x", matlab);
421  PrintVector(OType1y, My1, "OType1y_before_operation", matlab);
422  PrintVector(OType2x, Mx2, "OType2x", matlab);
423  PrintVector(OType2y, My2, "OType2y_before_operation", matlab);
424  }
425  TotalTestCount++;
426  OType1BLAS.AXPY(M1, OType1alpha, OType1x, incx1, OType1y, incy1);
427  OType2BLAS.AXPY(M2, OType2alpha, OType2x, incx2, OType2y, incy2);
428  if(debug)
429  {
430  PrintVector(OType1y, My1, "OType1y_after_operation", matlab);
431  PrintVector(OType2y, My2, "OType2y_after_operation", matlab);
432  }
433  if (CompareVectors(OType1y, My1, OType2y, My2, TOL)==0)
434  std::cout << "FAILED TEST!!!!!!" << std::endl;
435  GoodTestSubcount += CompareVectors(OType1y, My1, OType2y, My2, TOL);
436 
437  delete [] OType1x;
438  delete [] OType1y;
439  delete [] OType2x;
440  delete [] OType2y;
441  }
442  GoodTestCount += GoodTestSubcount;
443  if(verbose || debug) std::cout << "AXPY: " << GoodTestSubcount << " of " << AXPYTESTS << " tests were successful." << std::endl;
444  if(debug) std::cout << std::endl;
445  //--------------------------------------------------------------------------------
446  // End AXPY Tests
447  //--------------------------------------------------------------------------------
448 
449  //--------------------------------------------------------------------------------
450  // Begin COPY Tests
451  //--------------------------------------------------------------------------------
452  GoodTestSubcount = 0;
453  for(i = 0; i < COPYTESTS; i++)
454  {
455  incx1 = GetRandom(1, MVMAX);
456  incy1 = GetRandom(1, MVMAX);
457  incx2 = ConvertType( incx1, incx2 );
458  incy2 = ConvertType( incy1, incy2 );
459  M1 = GetRandom(MVMIN, MVMAX);
460  M2 = ConvertType( M1, M2 );
461  Mx1 = M1*std::abs(incx1);
462  My1 = M1*std::abs(incy1);
463  if (Mx1 == 0) { Mx1 = 1; }
464  if (My1 == 0) { My1 = 1; }
465  Mx2 = ConvertType( Mx1, Mx2 );
466  My2 = ConvertType( My1, My2 );
467  OType1x = new SType[Mx1];
468  OType1y = new SType[My1];
469  OType2x = new SType[Mx2];
470  OType2y = new SType[My2];
471  for(j1 = 0, j2 = 0; j1 < Mx1; j1++, j2++)
472  {
473  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
474  OType2x[j2] = OType1x[j1];
475  }
476  for(j1 = 0, j2 = 0; j1 < My1; j1++, j2++)
477  {
478  OType1y[j1] = GetRandom(-SCALARMAX, SCALARMAX);
479  OType2y[j2] = OType1y[j1];
480  }
481  if(debug)
482  {
483  std::cout << "Test #" << TotalTestCount << " -- COPY -- " << std::endl;
484  PrintVector(OType1x, Mx1, "OType1x", matlab);
485  PrintVector(OType1y, My1, "OType1y_before_operation", matlab);
486  PrintVector(OType2x, Mx2, "OType2x", matlab);
487  PrintVector(OType2y, My2, "OType2y_before_operation", matlab);
488  }
489  TotalTestCount++;
490  OType1BLAS.COPY(M1, OType1x, incx1, OType1y, incy1);
491  OType2BLAS.COPY(M2, OType2x, incx2, OType2y, incy2);
492  if(debug)
493  {
494  PrintVector(OType1y, My1, "OType1y_after_operation", matlab);
495  PrintVector(OType2y, My2, "OType2y_after_operation", matlab);
496  }
497  if (CompareVectors(OType1y, My1, OType2y, My2, TOL) == 0 )
498  std::cout << "FAILED TEST!!!!!!" << std::endl;
499  GoodTestSubcount += CompareVectors(OType1y, My1, OType2y, My2, TOL);
500 
501  delete [] OType1x;
502  delete [] OType1y;
503  delete [] OType2x;
504  delete [] OType2y;
505  }
506  GoodTestCount += GoodTestSubcount; if(verbose || debug) std::cout << "COPY: " << GoodTestSubcount << " of " << COPYTESTS << " tests were successful." << std::endl;
507  if(debug) std::cout << std::endl;
508  //--------------------------------------------------------------------------------
509  // End COPY Tests
510  //--------------------------------------------------------------------------------
511 
512  //--------------------------------------------------------------------------------
513  // Begin DOT Tests
514  //--------------------------------------------------------------------------------
515  GoodTestSubcount = 0;
516  for(i = 0; i < DOTTESTS; i++)
517  {
518  incx1 = GetRandom(1, MVMAX);
519  incy1 = GetRandom(1, MVMAX);
520  incx2 = ConvertType( incx1, incx2 );
521  incy2 = ConvertType( incy1, incy2 );
522  M1 = GetRandom(MVMIN, MVMAX);
523  M2 = ConvertType( M1, M2 );
524  Mx1 = M1*std::abs(incx1);
525  My1 = M1*std::abs(incy1);
526  if (Mx1 == 0) { Mx1 = 1; }
527  if (My1 == 0) { My1 = 1; }
528  Mx2 = ConvertType( Mx1, Mx2 );
529  My2 = ConvertType( My1, My2 );
530  OType1x = new SType[Mx1];
531  OType1y = new SType[My1];
532  OType2x = new SType[Mx2];
533  OType2y = new SType[My2];
534  for(j1 = 0, j2 = 0; j1 < Mx1; j1++, j2++)
535  {
536  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
537  OType2x[j2] = OType1x[j1];
538  }
539  for(j1 = 0, j2 = 0; j1 < My1; j1++, j2++)
540  {
541  OType1y[j1] = GetRandom(-SCALARMAX, SCALARMAX);
542  OType2y[j2] = OType1y[j1];
543  }
544  if(debug)
545  {
546  std::cout << "Test #" << TotalTestCount << " -- DOT -- " << std::endl;
547  PrintVector(OType1x, Mx1, "OType1x", matlab);
548  PrintVector(OType1y, My1, "OType1y", matlab);
549  PrintVector(OType2x, Mx2, "OType2x", matlab);
550  PrintVector(OType2y, My2, "OType2y", matlab);
551  }
552  TotalTestCount++;
553  OType1DOTresult = OType1BLAS.DOT(M1, OType1x, incx1, OType1y, incy1);
554  OType2DOTresult = OType2BLAS.DOT(M2, OType2x, incx2, OType2y, incy2);
555  if(debug)
556  {
557  std::cout << "OType1 DOT result: " << OType1DOTresult << std::endl;
558  std::cout << "OType2 DOT result: " << OType2DOTresult << std::endl;
559  }
560  if (CompareScalars(OType1DOTresult, OType2DOTresult, TOL) == 0) {
561  std::cout << "DOT test " << i+1 << " of " << DOTTESTS << " FAILED! "
562  << "SType = " << Teuchos::TypeNameTraits<SType>::name () << ". "
563  << "The two results are " << OType1DOTresult << " and "
564  << OType2DOTresult << ". incx1 = " << incx1 << ", incy1 = "
565  << incy1 << ", incx2 = " << incx2 << ", and incy2 = "
566  << incy2 << std::endl;
567  }
568 
569  GoodTestSubcount += CompareScalars(OType1DOTresult, OType2DOTresult, TOL);
570 
571  delete [] OType1x;
572  delete [] OType1y;
573  delete [] OType2x;
574  delete [] OType2y;
575  }
576  GoodTestCount += GoodTestSubcount;
577  if(verbose || debug) std::cout << "DOT: " << GoodTestSubcount << " of " << DOTTESTS << " tests were successful." << std::endl;
578  if(debug) std::cout << std::endl;
579  //--------------------------------------------------------------------------------
580  // End DOT Tests
581  //--------------------------------------------------------------------------------
582 
583  //--------------------------------------------------------------------------------
584  // Begin NRM2 Tests
585  //--------------------------------------------------------------------------------
586  GoodTestSubcount = 0;
587  for(i = 0; i < NRM2TESTS; i++)
588  {
589  incx1 = GetRandom(1, MVMAX);
590  incx2 = ConvertType( incx1, incx2 );
591  M1 = GetRandom(MVMIN, MVMAX);
592  M2 = ConvertType( M1, M2 );
593  OType1x = new SType[M1*incx1];
594  OType2x = new SType[M2*incx2];
595  for(j1 = 0, j2 = 0; j1 < M1*incx1; j1++, j2++)
596  {
597  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
598  OType2x[j2] = OType1x[j1];
599  }
600  if(debug)
601  {
602  std::cout << "Test #" << TotalTestCount << " -- NRM2 -- " << std::endl;
603  PrintVector(OType1x, M1*incx1, "OType1x", matlab);
604  PrintVector(OType2x, M2*incx2, "OType2x", matlab);
605  }
606  TotalTestCount++;
607  OType1NRM2result = OType1BLAS.NRM2(M1, OType1x, incx1);
608  OType2NRM2result = OType2BLAS.NRM2(M2, OType2x, incx2);
609  if(debug)
610  {
611  std::cout << "OType1 NRM2 result: " << OType1NRM2result << std::endl;
612  std::cout << "OType2 NRM2 result: " << OType2NRM2result << std::endl;
613  }
614  if (CompareScalars(OType1NRM2result, OType2NRM2result, TOL)==0)
615  std::cout << "FAILED TEST!!!!!!" << std::endl;
616  GoodTestSubcount += CompareScalars(OType1NRM2result, OType2NRM2result, TOL);
617 
618  delete [] OType1x;
619  delete [] OType2x;
620  }
621  GoodTestCount += GoodTestSubcount; if(verbose || debug) std::cout << "NRM2: " << GoodTestSubcount << " of " << NRM2TESTS << " tests were successful." << std::endl;
622  if(debug) std::cout << std::endl;
623  //--------------------------------------------------------------------------------
624  // End NRM2 Tests
625  //--------------------------------------------------------------------------------
626 
627  //--------------------------------------------------------------------------------
628  // Begin SCAL Tests
629  //--------------------------------------------------------------------------------
630  GoodTestSubcount = 0;
631  for(i = 0; i < SCALTESTS; i++)
632  {
633  // These will only test for the case that the increment is > 0, the
634  // templated case can handle when incx < 0, but the blas library doesn't
635  // seem to be able to on some machines.
636  incx1 = GetRandom(1, MVMAX);
637  incx2 = ConvertType( incx1, incx2 );
638  M1 = GetRandom(MVMIN, MVMAX);
639  M2 = ConvertType( M1, M2 );
640  OType1x = new SType[M1*incx1];
641  OType2x = new SType[M2*incx2];
642  for(j1 = 0, j2 = 0; j1 < M1*incx1; j1++, j2++)
643  {
644  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
645  OType2x[j2] = OType1x[j1];
646  }
647  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
648  OType2alpha = OType1alpha;
649  if(debug)
650  {
651  std::cout << "Test #" << TotalTestCount << " -- SCAL -- " << std::endl;
652  std::cout << "OType1alpha = " << OType1alpha << std::endl;
653  std::cout << "OType2alpha = " << OType2alpha << std::endl;
654  PrintVector(OType1x, M1*incx1, "OType1x_before_operation", matlab);
655  PrintVector(OType2x, M2*incx2, "OType2x_before_operation", matlab);
656  }
657  TotalTestCount++;
658  OType1BLAS.SCAL(M1, OType1alpha, OType1x, incx1);
659  OType2BLAS.SCAL(M2, OType2alpha, OType2x, incx2);
660  if(debug)
661  {
662  PrintVector(OType1x, M1*incx1, "OType1x_after_operation", matlab);
663  PrintVector(OType2x, M2*incx2, "OType2x_after_operation", matlab);
664  }
665  if (CompareVectors(OType1x, M1*incx1, OType2x, M2*incx2, TOL)==0)
666  std::cout << "FAILED TEST!!!!!!" << std::endl;
667  GoodTestSubcount += CompareVectors(OType1x, M1*incx1, OType2x, M2*incx2, TOL);
668 
669  delete [] OType1x;
670  delete [] OType2x;
671  }
672  GoodTestCount += GoodTestSubcount;
673  if(verbose || debug) std::cout << "SCAL: " << GoodTestSubcount << " of " << SCALTESTS << " tests were successful." << std::endl;
674  if(debug) std::cout << std::endl;
675  //--------------------------------------------------------------------------------
676  // End SCAL Tests
677  //--------------------------------------------------------------------------------
678 
679  //--------------------------------------------------------------------------------
680  // Begin IAMAX Tests
681  //--------------------------------------------------------------------------------
682  GoodTestSubcount = 0;
683  for(i = 0; i < IAMAXTESTS; i++)
684  {
685  incx1 = GetRandom(1, MVMAX);
686  incx2 = ConvertType( incx1, incx2 );
687  M1 = GetRandom(MVMIN, MVMAX);
688  M2 = ConvertType( M1, M2 );
689  OType1x = new SType[M1*incx1];
690  OType2x = new SType[M2*incx2];
691  for(j1 = 0, j2 = 0; j1 < M1*incx1; j1++, j2++)
692  {
693  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
694  OType2x[j2] = OType1x[j1];
695  }
696  if(debug)
697  {
698  std::cout << "Test #" << TotalTestCount << " -- IAMAX -- " << std::endl;
699  PrintVector(OType1x, M1*incx1, "OType1x", matlab);
700  PrintVector(OType2x, M2*incx2, "OType2x", matlab);
701  }
702  TotalTestCount++;
703  OType1IAMAXresult = OType1BLAS.IAMAX(M1, OType1x, incx1);
704  OType2IAMAXresult = OType2BLAS.IAMAX(M2, OType2x, incx2);
705  if(debug)
706  {
707  std::cout << "OType1 IAMAX result: " << OType1IAMAXresult << std::endl;
708  std::cout << "OType2 IAMAX result: " << OType2IAMAXresult << std::endl;
709  }
710  if (OType1IAMAXresult != OType2IAMAXresult)
711  std::cout << "FAILED TEST!!!!!!" << std::endl;
712  GoodTestSubcount += (OType1IAMAXresult == OType2IAMAXresult);
713 
714  delete [] OType1x;
715  delete [] OType2x;
716  }
717  GoodTestCount += GoodTestSubcount;
718  if(verbose || debug) std::cout << "IAMAX: " << GoodTestSubcount << " of " << IAMAXTESTS << " tests were successful." << std::endl;
719  if(debug) std::cout << std::endl;
720  //--------------------------------------------------------------------------------
721  // End IAMAX Tests
722  //--------------------------------------------------------------------------------
723 
724  //--------------------------------------------------------------------------------
725  // BEGIN LEVEL 2 BLAS TESTS
726  //--------------------------------------------------------------------------------
727  // Begin GEMV Tests
728  //--------------------------------------------------------------------------------
729  GoodTestSubcount = 0;
730  for(i = 0; i < GEMVTESTS; i++)
731  {
732  // The parameters used to construct the test problem are chosen to be
733  // valid parameters, so the GEMV routine won't bomb out.
734  incx1 = GetRandom(1, MVMAX);
735  while (incx1 == 0) {
736  incx1 = GetRandom(1, MVMAX);
737  }
738  incy1 = GetRandom(1, MVMAX);
739  while (incy1 == 0) {
740  incy1 = GetRandom(1, MVMAX);
741  }
742  incx2 = ConvertType( incx1, incx2 );
743  incy2 = ConvertType( incy1, incy2 );
744  M1 = GetRandom(MVMIN, MVMAX);
745  N1 = GetRandom(MVMIN, MVMAX);
746  M2 = ConvertType( M1, M2 );
747  N2 = ConvertType( N1, N2 );
748 
749  TRANS = RandomTRANS();
750  OType1 M2_1 = 0, N2_1 = 0;
751  if (Teuchos::ETranspChar[TRANS] == 'N') {
752  M2_1 = M1*std::abs(incy1);
753  N2_1 = N1*std::abs(incx1);
754  } else {
755  M2_1 = N1*std::abs(incy1);
756  N2_1 = M1*std::abs(incx1);
757  }
758  OType2 M2_2 = 0;
759  OType2 N2_2 = 0;
760  M2_2 = ConvertType( M2_1, M2_2 );
761  N2_2 = ConvertType( N2_1, N2_2 );
762 
763  LDA1 = GetRandom(MVMIN, MVMAX);
764  while (LDA1 < M1) {
765  LDA1 = GetRandom(MVMIN, MVMAX);
766  }
767  LDA2 = ConvertType( LDA1, LDA2 );
768 
769  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
770  OType1beta = GetRandom(-SCALARMAX, SCALARMAX);
771  OType2alpha = OType1alpha;
772  OType2beta = OType1beta;
773 
774  OType1A = new SType[LDA1 * N1];
775  OType1x = new SType[N2_1];
776  OType1y = new SType[M2_1];
777  OType2A = new SType[LDA2 * N2];
778  OType2x = new SType[N2_2];
779  OType2y = new SType[M2_2];
780 
781  for(j1 = 0, j2 = 0; j1 < LDA1 * N1; j1++, j2++)
782  {
783  OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX);
784  OType2A[j2] = OType1A[j1];
785  }
786  for(j1 = 0, j2 = 0; j1 < N2_1; j1++, j2++)
787  {
788  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
789  OType2x[j2] = OType1x[j1];
790  }
791  for(j1 = 0, j2 = 0; j1 < M2_1; j1++, j2++)
792  {
793  OType1y[j1] = GetRandom(-SCALARMAX, SCALARMAX);
794  OType2y[j2] = OType1y[j1];
795  }
796  if(debug)
797  {
798  std::cout << "Test #" << TotalTestCount << " -- GEMV -- " << std::endl;
799  std::cout << "TRANS = " << Teuchos::ETranspChar[TRANS] << std::endl;
800  std::cout << "OType1alpha = " << OType1alpha << std::endl;
801  std::cout << "OType2alpha = " << OType2alpha << std::endl;
802  std::cout << "OType1beta = " << OType1beta << std::endl;
803  std::cout << "OType2beta = " << OType2beta << std::endl;
804  PrintMatrix(OType1A, M1, N1, LDA1, "OType1A", matlab);
805  PrintVector(OType1x, N2_1, "OType1x", matlab);
806  PrintVector(OType1y, M2_1, "OType1y_before_operation", matlab);
807  PrintMatrix(OType2A, M2, N2, LDA2, "OType2A", matlab);
808  PrintVector(OType2x, N2_2, "OType2x", matlab);
809  PrintVector(OType2y, M2_2, "OType2y_before_operation", matlab);
810  }
811  TotalTestCount++;
812  OType1BLAS.GEMV(TRANS, M1, N1, OType1alpha, OType1A, LDA1, OType1x, incx1, OType1beta, OType1y, incy1);
813  OType2BLAS.GEMV(TRANS, M2, N2, OType2alpha, OType2A, LDA2, OType2x, incx2, OType2beta, OType2y, incy2);
814  if(debug)
815  {
816  PrintVector(OType1y, M2_1, "OType1y_after_operation", matlab);
817  PrintVector(OType2y, M2_2, "OType2y_after_operation", matlab);
818  }
819  if (CompareVectors(OType1y, M2_1, OType2y, M2_2, TOL)==0)
820  std::cout << "FAILED TEST!!!!!!" << std::endl;
821  GoodTestSubcount += CompareVectors(OType1y, M2_1, OType2y, M2_2, TOL);
822 
823  delete [] OType1A;
824  delete [] OType1x;
825  delete [] OType1y;
826  delete [] OType2A;
827  delete [] OType2x;
828  delete [] OType2y;
829  }
830  GoodTestCount += GoodTestSubcount;
831  if(verbose || debug) std::cout << "GEMV: " << GoodTestSubcount << " of " << GEMVTESTS << " tests were successful." << std::endl;
832  if(debug) std::cout << std::endl;
833  //--------------------------------------------------------------------------------
834  // End GEMV Tests
835  //--------------------------------------------------------------------------------
836 
837  //--------------------------------------------------------------------------------
838  // Begin TRMV Tests
839  //--------------------------------------------------------------------------------
840  GoodTestSubcount = 0;
841  for(i = 0; i < TRMVTESTS; i++)
842  {
843  UPLO = RandomUPLO();
844  TRANSA = RandomTRANS();
845 
846  // Since the entries are integers, we don't want to use the unit diagonal feature,
847  // this creates ill-conditioned, nearly-singular matrices.
848  //DIAG = RandomDIAG();
849  DIAG = Teuchos::NON_UNIT_DIAG;
850 
851  N1 = GetRandom(MVMIN, MVMAX);
852  N2 = ConvertType( N1, N2 );
853  incx1 = GetRandom(1, MVMAX);
854  while (incx1 == 0) {
855  incx1 = GetRandom(1, MVMAX);
856  }
857  incx2 = ConvertType( incx1, incx2 );
858  OType1x = new SType[N1*std::abs(incx1)];
859  OType2x = new SType[N2*std::abs(incx2)];
860 
861  for(j1 = 0, j2 = 0; j1 < N1*std::abs(incx1); j1++, j2++)
862  {
863  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
864  OType2x[j2] = OType1x[j1];
865  }
866 
867  LDA1 = GetRandom(MVMIN, MVMAX);
868  while (LDA1 < N1) {
869  LDA1 = GetRandom(MVMIN, MVMAX);
870  }
871  LDA2 = ConvertType( LDA1, LDA2 );
872  OType1A = new SType[LDA1 * N1];
873  OType2A = new SType[LDA2 * N2];
874 
875  for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++)
876  {
877  if(Teuchos::EUploChar[UPLO] == 'U') {
878  // The operator is upper triangular, make sure that the entries are
879  // only in the upper triangular part of A and the diagonal is non-zero.
880  for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++)
881  {
882  if(k1 < j1) {
883  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
884  } else {
885  OType1A[j1*LDA1+k1] = STypezero;
886  }
887  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
888  if(k1 == j1) {
889  if (Teuchos::EDiagChar[DIAG] == 'N') {
890  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
891  while (OType1A[j1*LDA1+k1] == STypezero) {
892  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
893  }
894  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
895  } else {
896  OType1A[j1*LDA1+k1] = STypeone;
897  OType2A[j2*LDA2+k2] = STypeone;
898  }
899  }
900  }
901  } else {
902  // The operator is lower triangular, make sure that the entries are
903  // only in the lower triangular part of A and the diagonal is non-zero.
904  for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++)
905  {
906  if(k1 > j1) {
907  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
908  } else {
909  OType1A[j1*LDA1+k1] = STypezero;
910  }
911  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
912  if(k1 == j1) {
913  if (Teuchos::EDiagChar[DIAG] == 'N') {
914  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
915  while (OType1A[j1*LDA1+k1] == STypezero) {
916  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
917  }
918  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
919  } else {
920  OType1A[j1*LDA1+k1] = STypeone;
921  OType2A[j2*LDA2+k2] = STypeone;
922  }
923  }
924  } // end for(k=0 ...
925  } // end if(UPLO == 'U') ...
926  } // end for(j=0 ... for(j = 0; j < N*N; j++)
927 
928  if(debug)
929  {
930  std::cout << "Test #" << TotalTestCount << " -- TRMV -- " << std::endl;
931  std::cout << "UPLO = " << Teuchos::EUploChar[UPLO] << "\t"
932  << "TRANSA = " << Teuchos::ETranspChar[TRANSA] << "\t"
933  << "DIAG = " << Teuchos::EDiagChar[DIAG] << std::endl;
934  PrintMatrix(OType1A, N1, N1, LDA1,"OType1A", matlab);
935  PrintVector(OType1x, N1*incx1, "OType1x_before_operation", matlab);
936  PrintMatrix(OType2A, N2, N2, LDA2, "OType2A", matlab);
937  PrintVector(OType2x, N2*incx2, "OType2x_before_operation", matlab);
938  }
939  TotalTestCount++;
940  OType1BLAS.TRMV(UPLO, TRANSA, DIAG, N1, OType1A, LDA1, OType1x, incx1);
941  OType2BLAS.TRMV(UPLO, TRANSA, DIAG, N2, OType2A, LDA2, OType2x, incx2);
942  if(debug)
943  {
944  PrintVector(OType1x, N1*incx1, "OType1x_after_operation", matlab);
945  PrintVector(OType2x, N2*incx2, "OType2x_after_operation", matlab);
946  }
947  if (CompareVectors(OType1x, std::abs(N1*incx1), OType2x, std::abs(N2*incx2), TOL)==0)
948  std::cout << "FAILED TEST!!!!!!" << std::endl;
949  GoodTestSubcount += CompareVectors(OType1x, std::abs(N1*incx1), OType2x, std::abs(N2*incx2), TOL);
950 
951  delete [] OType1A;
952  delete [] OType1x;
953  delete [] OType2A;
954  delete [] OType2x;
955  }
956  GoodTestCount += GoodTestSubcount;
957  if(verbose || debug) std::cout << "TRMV: " << GoodTestSubcount << " of " << TRMVTESTS << " tests were successful." << std::endl;
958  if(debug) std::cout << std::endl;
959  //--------------------------------------------------------------------------------
960  // End TRMV Tests
961  //--------------------------------------------------------------------------------
962 
963  //--------------------------------------------------------------------------------
964  // Begin GER Tests
965  //--------------------------------------------------------------------------------
966  GoodTestSubcount = 0;
967  for(i = 0; i < GERTESTS; i++)
968  {
969  incx1 = GetRandom(1, MVMAX);
970  while (incx1 == 0) {
971  incx1 = GetRandom(1, MVMAX);
972  }
973  incy1 = GetRandom(1, MVMAX);
974  while (incy1 == 0) {
975  incy1 = GetRandom(1, MVMAX);
976  }
977  incx2 = ConvertType( incx1, incx2 );
978  incy2 = ConvertType( incy1, incy2 );
979  M1 = GetRandom(MVMIN, MVMAX);
980  N1 = GetRandom(MVMIN, MVMAX);
981  M2 = ConvertType( M1, M2 );
982  N2 = ConvertType( N1, N2 );
983 
984  LDA1 = GetRandom(MVMIN, MVMAX);
985  while (LDA1 < M1) {
986  LDA1 = GetRandom(MVMIN, MVMAX);
987  }
988  LDA2 = ConvertType( LDA1, LDA2 );
989 
990  OType1A = new SType[LDA1 * N1];
991  OType1x = new SType[M1*std::abs(incx1)];
992  OType1y = new SType[N1*std::abs(incy1)];
993  OType2A = new SType[LDA2 * N2];
994  OType2x = new SType[M2*std::abs(incx2)];
995  OType2y = new SType[N2*std::abs(incy2)];
996  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
997  OType2alpha = OType1alpha;
998  for(j1 = 0, j2 = 0; j1 < LDA1 * N1; j1++, j2++)
999  {
1000  OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1001  OType2A[j2] = OType1A[j1];
1002  }
1003  for(j1 = 0, j2 = 0; j1 < std::abs(M1*incx1); j1++, j2++)
1004  {
1005  OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1006  OType2x[j2] = OType1x[j1];
1007  }
1008  for(j1 = 0, j2 = 0; j1 < std::abs(N1*incy1); j1++, j2++)
1009  {
1010  OType1y[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1011  OType2y[j2] = OType1y[j1];
1012  }
1013  if(debug)
1014  {
1015  std::cout << "Test #" << TotalTestCount << " -- GER -- " << std::endl;
1016  std::cout << "OType1alpha = " << OType1alpha << std::endl;
1017  std::cout << "OType2alpha = " << OType2alpha << std::endl;
1018  PrintMatrix(OType1A, M1, N1, LDA1,"OType1A_before_operation", matlab);
1019  PrintVector(OType1x, std::abs(M1*incx1), "OType1x", matlab);
1020  PrintVector(OType1y, std::abs(N1*incy1), "OType1y", matlab);
1021  PrintMatrix(OType2A, M2, N2, LDA2,"OType2A_before_operation", matlab);
1022  PrintVector(OType2x, std::abs(M2*incx2), "OType2x", matlab);
1023  PrintVector(OType2y, std::abs(N2*incy2), "OType2y", matlab);
1024  }
1025  TotalTestCount++;
1026  OType1BLAS.GER(M1, N1, OType1alpha, OType1x, incx1, OType1y, incy1, OType1A, LDA1);
1027  OType2BLAS.GER(M2, N2, OType2alpha, OType2x, incx2, OType2y, incy2, OType2A, LDA2);
1028  if(debug)
1029  {
1030  PrintMatrix(OType1A, M1, N1, LDA1, "OType1A_after_operation", matlab);
1031  PrintMatrix(OType2A, M2, N2, LDA2, "OType2A_after_operation", matlab);
1032  }
1033  if (CompareMatrices(OType1A, M1, N1, LDA1, OType2A, M2, N2, LDA2, TOL)==0)
1034  std::cout << "FAILED TEST!!!!!!" << std::endl;
1035  GoodTestSubcount += CompareMatrices(OType1A, M1, N1, LDA1, OType2A, M2, N2, LDA2, TOL);
1036 
1037  delete [] OType1A;
1038  delete [] OType1x;
1039  delete [] OType1y;
1040  delete [] OType2A;
1041  delete [] OType2x;
1042  delete [] OType2y;
1043  }
1044  GoodTestCount += GoodTestSubcount;
1045  if(verbose || debug) std::cout << "GER: " << GoodTestSubcount << " of " << GERTESTS << " tests were successful." << std::endl;
1046  if(debug) std::cout << std::endl;
1047  //--------------------------------------------------------------------------------
1048  // End GER Tests
1049  //--------------------------------------------------------------------------------
1050 
1051  //--------------------------------------------------------------------------------
1052  // BEGIN LEVEL 3 BLAS TESTS
1053  //--------------------------------------------------------------------------------
1054  // Begin GEMM Tests
1055  //--------------------------------------------------------------------------------
1056  GoodTestSubcount = 0;
1057  for(i = 0; i < GEMMTESTS; i++)
1058  {
1059  TRANSA = RandomTRANS();
1060  TRANSB = RandomTRANS();
1061  M1 = GetRandom(MVMIN, MVMAX);
1062  N1 = GetRandom(MVMIN, MVMAX);
1063  P1 = GetRandom(MVMIN, MVMAX);
1064  M2 = ConvertType( M1, M2 );
1065  N2 = ConvertType( N1, N2 );
1066  P2 = ConvertType( P1, P2 );
1067 
1068  if(debug) {
1069  std::cout << "Test #" << TotalTestCount << " -- GEMM -- " << std::endl;
1070  std::cout << "TRANSA = " << Teuchos::ETranspChar[TRANSA] << "\t"
1071  << "TRANSB = " << Teuchos::ETranspChar[TRANSB] << std::endl;
1072  }
1073  LDA1 = GetRandom(MVMIN, MVMAX);
1074  if (Teuchos::ETranspChar[TRANSA] == 'N') {
1075  while (LDA1 < M1) { LDA1 = GetRandom(MVMIN, MVMAX); }
1076  LDA2 = ConvertType( LDA1, LDA2 );
1077  OType1A = new SType[LDA1 * P1];
1078  OType2A = new SType[LDA2 * P2];
1079  for(j1 = 0, j2 = 0; j1 < LDA1 * P1; j1++, j2++)
1080  {
1081  OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1082  OType2A[j2] = OType1A[j1];
1083  }
1084  if (debug) {
1085  PrintMatrix(OType1A, M1, P1, LDA1, "OType1A", matlab);
1086  PrintMatrix(OType2A, M2, P2, LDA2, "OType2A", matlab);
1087  }
1088  } else {
1089  while (LDA1 < P1) { LDA1 = GetRandom(MVMIN, MVMAX); }
1090  LDA2 = ConvertType( LDA1, LDA2 );
1091  OType1A = new SType[LDA1 * M1];
1092  OType2A = new SType[LDA1 * M1];
1093  for(j1 = 0, j2 = 0; j1 < LDA1 * M1; j1++, j2++)
1094  {
1095  OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1096  OType2A[j2] = OType1A[j1];
1097  }
1098  if (debug) {
1099  PrintMatrix(OType1A, P1, M1, LDA1, "OType1A", matlab);
1100  PrintMatrix(OType2A, P2, M2, LDA2, "OType2A", matlab);
1101  }
1102  }
1103 
1104  LDB1 = GetRandom(MVMIN, MVMAX);
1105  if (Teuchos::ETranspChar[TRANSB] == 'N') {
1106  while (LDB1 < P1) { LDB1 = GetRandom(MVMIN, MVMAX); }
1107  LDB2 = ConvertType( LDB1, LDB2 );
1108  OType1B = new SType[LDB1 * N1];
1109  OType2B = new SType[LDB2 * N2];
1110  for(j1 = 0, j2 = 0; j1 < LDB1 * N1; j1++, j2++)
1111  {
1112  OType1B[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1113  OType2B[j2] = OType1B[j1];
1114  }
1115  if (debug) {
1116  PrintMatrix(OType1B, P1, N1, LDB1,"OType1B", matlab);
1117  PrintMatrix(OType2B, P2, N2, LDB2,"OType2B", matlab);
1118  }
1119  } else {
1120  while (LDB1 < N1) { LDB1 = GetRandom(MVMIN, MVMAX); }
1121  LDB2 = ConvertType( LDB1, LDB2 );
1122  OType1B = new SType[LDB1 * P1];
1123  OType2B = new SType[LDB2 * P2];
1124  for(j1 = 0, j2 = 0; j1 < LDB1 * P1; j1++, j2++)
1125  {
1126  OType1B[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1127  OType2B[j2] = OType1B[j1];
1128  }
1129  if (debug) {
1130  PrintMatrix(OType1B, N1, P1, LDB1,"OType1B", matlab);
1131  PrintMatrix(OType2B, N2, P2, LDB2,"OType2B", matlab);
1132  }
1133  }
1134 
1135  LDC1 = GetRandom(MVMIN, MVMAX);
1136  while (LDC1 < M1) { LDC1 = GetRandom(MVMIN, MVMAX); }
1137  LDC2 = ConvertType( LDC1, LDC2 );
1138  OType1C = new SType[LDC1 * N1];
1139  OType2C = new SType[LDC2 * N2];
1140  for(j1 = 0, j2 = 0; j1 < LDC1 * N1; j1++, j2++) {
1141  OType1C[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1142  OType2C[j2] = OType1C[j1];
1143  }
1144  if(debug)
1145  {
1146  PrintMatrix(OType1C, M1, N1, LDC1, "OType1C_before_operation", matlab);
1147  PrintMatrix(OType2C, M2, N2, LDC2, "OType2C_before_operation", matlab);
1148  }
1149 
1150  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
1151  OType1beta = GetRandom(-SCALARMAX, SCALARMAX);
1152  OType2alpha = OType1alpha;
1153  OType2beta = OType1beta;
1154 
1155  TotalTestCount++;
1156  OType1BLAS.GEMM(TRANSA, TRANSB, M1, N1, P1, OType1alpha, OType1A, LDA1, OType1B, LDB1, OType1beta, OType1C, LDC1);
1157  OType2BLAS.GEMM(TRANSA, TRANSB, M2, N2, P2, OType2alpha, OType2A, LDA2, OType2B, LDB2, OType2beta, OType2C, LDC2);
1158  if(debug)
1159  {
1160  std::cout << "M1="<<M1 << "\t" << "N1="<<N1 << "\t" << "P1 = " << P1
1161  << "\t" << "LDA1="<<LDA1 << "\t" << "LDB1="<<LDB1 << "\t" << "LDC1=" << LDC1 << std::endl;
1162  std::cout << "M2="<<M2 << "\t" << "N2="<<N2 << "\t" << "P2 = " << P2
1163  << "\t" << "LDA2="<<LDA2 << "\t" << "LDB2="<<LDB2 << "\t" << "LDC2=" << LDC2 << std::endl;
1164  std::cout << "OType1alpha = " << OType1alpha << std::endl;
1165  std::cout << "OType2alpha = " << OType2alpha << std::endl;
1166  std::cout << "OType1beta = " << OType1beta << std::endl;
1167  std::cout << "OType2beta = " << OType2beta << std::endl;
1168  PrintMatrix(OType1C, M1, N1, LDC1, "OType1C_after_operation", matlab);
1169  PrintMatrix(OType2C, M2, N2, LDC2, "OType2C_after_operation", matlab);
1170  }
1171  if (CompareMatrices(OType1C, M1, N1, LDC1, OType2C, M2, N2, LDC2, TOL)==0)
1172  std::cout << "FAILED TEST!!!!!!" << std::endl;
1173  GoodTestSubcount += CompareMatrices(OType1C, M1, N1, LDC1, OType2C, M2, N2, LDC2, TOL);
1174 
1175  delete [] OType1A;
1176  delete [] OType1B;
1177  delete [] OType1C;
1178  delete [] OType2A;
1179  delete [] OType2B;
1180  delete [] OType2C;
1181  }
1182  GoodTestCount += GoodTestSubcount;
1183  if(verbose || debug) std::cout << "GEMM: " << GoodTestSubcount << " of " << GEMMTESTS << " tests were successful." << std::endl;
1184  if(debug) std::cout << std::endl;
1185  //--------------------------------------------------------------------------------
1186  // End GEMM Tests
1187  //--------------------------------------------------------------------------------
1188 
1189  //--------------------------------------------------------------------------------
1190  // Begin SYMM Tests
1191  //--------------------------------------------------------------------------------
1192  GoodTestSubcount = 0;
1193  for(i = 0; i < SYMMTESTS; i++)
1194  {
1195  M1 = GetRandom(MVMIN, MVMAX);
1196  N1 = GetRandom(MVMIN, MVMAX);
1197  M2 = ConvertType( M1, M2 );
1198  N2 = ConvertType( N1, N2 );
1199  SIDE = RandomSIDE();
1200  UPLO = RandomUPLO();
1201 
1202  LDA1 = GetRandom(MVMIN, MVMAX);
1203  if(Teuchos::ESideChar[SIDE] == 'L') {
1204  while (LDA1 < M1) { LDA1 = GetRandom(MVMIN, MVMAX); }
1205  LDA2 = ConvertType( LDA1, LDA2 );
1206  OType1A = new SType[LDA1 * M1];
1207  OType2A = new SType[LDA2 * M2];
1208  for(j1 = 0, j2 = 0; j1 < LDA1 * M1; j1++, j2++) {
1209  OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1210  OType2A[j2] = OType1A[j1];
1211  }
1212  } else {
1213  while (LDA1 < N1) { LDA1 = GetRandom(MVMIN, MVMAX); }
1214  LDA2 = ConvertType( LDA1, LDA2 );
1215  OType1A = new SType[LDA1 * N1];
1216  OType2A = new SType[LDA2 * N2];
1217  for(j1 = 0, j2 = 0; j1 < LDA1 * N1; j1++, j2++) {
1218  OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1219  OType2A[j2] = OType1A[j1];
1220  }
1221  }
1222 
1223  LDB1 = GetRandom(MVMIN, MVMAX);
1224  while (LDB1 < M1) { LDB1 = GetRandom(MVMIN, MVMAX); }
1225  LDB2 = ConvertType( LDB1, LDB2 );
1226  OType1B = new SType[LDB1 * N1];
1227  OType2B = new SType[LDB2 * N2];
1228  for(j1 = 0, j2 = 0; j1 < LDB1 * N1; j1++, j2++) {
1229  OType1B[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1230  OType2B[j2] = OType1B[j1];
1231  }
1232 
1233  LDC1 = GetRandom(MVMIN, MVMAX);
1234  while (LDC1 < M1) { LDC1 = GetRandom(MVMIN, MVMAX); }
1235  LDC2 = ConvertType( LDC1, LDC2 );
1236  OType1C = new SType[LDC1 * N1];
1237  OType2C = new SType[LDC2 * N2];
1238  for(j1 = 0, j2 = 0; j1 < LDC1 * N1; j1++, j2++) {
1239  OType1C[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1240  OType2C[j2] = OType1C[j1];
1241  }
1242 
1243  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
1244  OType1beta = GetRandom(-SCALARMAX, SCALARMAX);
1245  OType2alpha = OType1alpha;
1246  OType2beta = OType1beta;
1247  if(debug)
1248  {
1249  std::cout << "Test #" << TotalTestCount << " -- SYMM -- " << std::endl;
1250  std::cout << "SIDE = " << Teuchos::ESideChar[SIDE] << "\t"
1251  << "UPLO = " << Teuchos::EUploChar[UPLO] << std::endl;
1252  std::cout << "OType1alpha = " << OType1alpha << std::endl;
1253  std::cout << "OType2alpha = " << OType2alpha << std::endl;
1254  std::cout << "OType1beta = " << OType1beta << std::endl;
1255  std::cout << "OType2beta = " << OType2beta << std::endl;
1256  if (Teuchos::ESideChar[SIDE] == 'L') {
1257  PrintMatrix(OType1A, M1, M1, LDA1,"OType1A", matlab);
1258  PrintMatrix(OType2A, M2, M2, LDA2,"OType2A", matlab);
1259  } else {
1260  PrintMatrix(OType1A, N1, N1, LDA1, "OType1A", matlab);
1261  PrintMatrix(OType2A, N2, N2, LDA2, "OType2A", matlab);
1262  }
1263  PrintMatrix(OType1B, M1, N1, LDB1,"OType1B", matlab);
1264  PrintMatrix(OType1C, M1, N1, LDC1,"OType1C_before_operation", matlab);
1265  PrintMatrix(OType2B, M2, N2, LDB2,"OType2B", matlab);
1266  PrintMatrix(OType2C, M2, N2, LDC2,"OType2C_before_operation", matlab);
1267  }
1268  TotalTestCount++;
1269 
1270  OType1BLAS.SYMM(SIDE, UPLO, M1, N1, OType1alpha, OType1A, LDA1, OType1B, LDB1, OType1beta, OType1C, LDC1);
1271  OType2BLAS.SYMM(SIDE, UPLO, M2, N2, OType2alpha, OType2A, LDA2, OType2B, LDB2, OType2beta, OType2C, LDC2);
1272  if(debug)
1273  {
1274  PrintMatrix(OType1C, M1, N1, LDC1,"OType1C_after_operation", matlab);
1275  PrintMatrix(OType2C, M2, N2, LDC2,"OType2C_after_operation", matlab);
1276  }
1277  if (CompareMatrices(OType1C, M1, N1, LDC1, OType2C, M2, N2, LDC2, TOL)==0)
1278  std::cout << "FAILED TEST!!!!!!" << std::endl;
1279  GoodTestSubcount += CompareMatrices(OType1C, M1, N1, LDC1, OType2C, M2, N2, LDC2, TOL);
1280 
1281  delete [] OType1A;
1282  delete [] OType1B;
1283  delete [] OType1C;
1284  delete [] OType2A;
1285  delete [] OType2B;
1286  delete [] OType2C;
1287  }
1288  GoodTestCount += GoodTestSubcount;
1289  if(verbose || debug) std::cout << "SYMM: " << GoodTestSubcount << " of " << SYMMTESTS << " tests were successful." << std::endl;
1290  if(debug) std::cout << std::endl;
1291  //--------------------------------------------------------------------------------
1292  // End SYMM Tests
1293  //--------------------------------------------------------------------------------
1294 
1295  //--------------------------------------------------------------------------------
1296  // Begin SYRK Tests
1297  //--------------------------------------------------------------------------------
1298  GoodTestSubcount = 0;
1299  for(i = 0; i < SYRKTESTS; i++)
1300  {
1301  N1 = GetRandom(MVMIN, MVMAX);
1302  K1 = GetRandom(MVMIN, MVMAX);
1303  while (K1 > N1) { K1 = GetRandom(MVMIN, MVMAX); }
1304  N2 = ConvertType( N1, N2 );
1305  K2 = ConvertType( K1, K2 );
1306 
1307  UPLO = RandomUPLO();
1308  TRANS = RandomTRANS();
1309 #ifdef HAVE_TEUCHOS_COMPLEX
1310  while (TRANS == Teuchos::CONJ_TRANS) { TRANS = RandomTRANS(); }
1311 #endif
1312 
1313  LDA1 = GetRandom(MVMIN, MVMAX);
1314  if(Teuchos::ETranspChar[TRANS] == 'N') {
1315  while (LDA1 < N1) { LDA1 = GetRandom(MVMIN, MVMAX); }
1316  LDA2 = ConvertType( LDA1, LDA2 );
1317  OType1A = new SType[LDA1 * K1];
1318  OType2A = new SType[LDA2 * K2];
1319  for(j1 = 0, j2 = 0; j1 < LDA1 * K1; j1++, j2++) {
1320  OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1321  OType2A[j2] = OType1A[j1];
1322  }
1323  } else {
1324  while (LDA1 < K1) { LDA1 = GetRandom(MVMIN, MVMAX); }
1325  LDA2 = ConvertType( LDA1, LDA2 );
1326  OType1A = new SType[LDA1 * N1];
1327  OType2A = new SType[LDA2 * N2];
1328  for(j1 = 0, j2 = 0; j1 < LDA1 * N1; j1++, j2++) {
1329  OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX);
1330  OType2A[j2] = OType1A[j1];
1331  }
1332  }
1333 
1334  LDC1 = GetRandom(MVMIN, MVMAX);
1335  while (LDC1 < N1) { LDC1 = GetRandom(MVMIN, MVMAX); }
1336  LDC2 = ConvertType( LDC1, LDC2 );
1337  OType1C = new SType[LDC1 * N1];
1338  OType2C = new SType[LDC2 * N2];
1339  for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++) {
1340 
1341  if(Teuchos::EUploChar[UPLO] == 'U') {
1342  // The operator is upper triangular, make sure that the entries are
1343  // only in the upper triangular part of C.
1344  for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++)
1345  {
1346  if(k1 <= j1) {
1347  OType1C[j1*LDC1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1348  } else {
1349  OType1C[j1*LDC1+k1] = STypezero;
1350  }
1351  OType2C[j2*LDC2+k2] = OType1C[j1*LDC1+k1];
1352  }
1353  }
1354  else {
1355  for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++)
1356  {
1357  if(k1 >= j1) {
1358  OType1C[j1*LDC1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1359  } else {
1360  OType1C[j1*LDC1+k1] = STypezero;
1361  }
1362  OType2C[j2*LDC2+k2] = OType1C[j1*LDC1+k1];
1363  }
1364  }
1365  }
1366 
1367  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
1368  OType1beta = GetRandom(-SCALARMAX, SCALARMAX);
1369  OType2alpha = OType1alpha;
1370  OType2beta = OType1beta;
1371  if(debug)
1372  {
1373  std::cout << "Test #" << TotalTestCount << " -- SYRK -- " << std::endl;
1374  std::cout << "UPLO = " << Teuchos::EUploChar[UPLO] << "\t"
1375  << "TRANS = " << Teuchos::ETranspChar[TRANS] << std::endl;
1376  std::cout << "N1="<<N1 << "\t" << "K1 = " << K1
1377  << "\t" << "LDA1="<<LDA1 << "\t" << "LDC1=" << LDC1 << std::endl;
1378  std::cout << "N2="<<N2 << "\t" << "K2 = " << K2
1379  << "\t" << "LDA2="<<LDA2 << "\t" << "LDC2=" << LDC2 << std::endl;
1380  std::cout << "OType1alpha = " << OType1alpha << std::endl;
1381  std::cout << "OType2alpha = " << OType2alpha << std::endl;
1382  std::cout << "OType1beta = " << OType1beta << std::endl;
1383  std::cout << "OType2beta = " << OType2beta << std::endl;
1384  if (Teuchos::ETranspChar[TRANS] == 'N') {
1385  PrintMatrix(OType1A, N1, K1, LDA1,"OType1A", matlab);
1386  PrintMatrix(OType2A, N2, K2, LDA2,"OType2A", matlab);
1387  } else {
1388  PrintMatrix(OType1A, K1, N1, LDA1, "OType1A", matlab);
1389  PrintMatrix(OType2A, K2, N2, LDA2, "OType2A", matlab);
1390  }
1391  PrintMatrix(OType1C, N1, N1, LDC1,"OType1C_before_operation", matlab);
1392  PrintMatrix(OType2C, N2, N2, LDC2,"OType2C_before_operation", matlab);
1393  }
1394  TotalTestCount++;
1395 
1396  OType1BLAS.SYRK(UPLO, TRANS, N1, K1, OType1alpha, OType1A, LDA1, OType1beta, OType1C, LDC1);
1397  OType2BLAS.SYRK(UPLO, TRANS, N2, K2, OType2alpha, OType2A, LDA2, OType2beta, OType2C, LDC2);
1398  if(debug)
1399  {
1400  PrintMatrix(OType1C, N1, N1, LDC1,"OType1C_after_operation", matlab);
1401  PrintMatrix(OType2C, N2, N2, LDC2,"OType2C_after_operation", matlab);
1402  }
1403  if (CompareMatrices(OType1C, N1, N1, LDC1, OType2C, N2, N2, LDC2, TOL)==0)
1404  std::cout << "FAILED TEST!!!!!!" << std::endl;
1405  GoodTestSubcount += CompareMatrices(OType1C, N1, N1, LDC1, OType2C, N2, N2, LDC2, TOL);
1406 
1407  delete [] OType1A;
1408  delete [] OType1C;
1409  delete [] OType2A;
1410  delete [] OType2C;
1411  }
1412  GoodTestCount += GoodTestSubcount;
1413  if(verbose || debug) std::cout << "SYRK: " << GoodTestSubcount << " of " << SYRKTESTS << " tests were successful." << std::endl;
1414  if(debug) std::cout << std::endl;
1415  //--------------------------------------------------------------------------------
1416  // End SYRK Tests
1417  //--------------------------------------------------------------------------------
1418 
1419  //--------------------------------------------------------------------------------
1420  // Begin TRMM Tests
1421  //--------------------------------------------------------------------------------
1422  GoodTestSubcount = 0;
1423  for(i = 0; i < TRMMTESTS; i++)
1424  {
1425  M1 = GetRandom(MVMIN, MVMAX);
1426  N1 = GetRandom(MVMIN, MVMAX);
1427  M2 = ConvertType( M1, M2 );
1428  N2 = ConvertType( N1, N2 );
1429 
1430  LDB1 = GetRandom(MVMIN, MVMAX);
1431  while (LDB1 < M1) {
1432  LDB1 = GetRandom(MVMIN, MVMAX);
1433  }
1434  LDB2 = ConvertType( LDB1, LDB2 );
1435 
1436  OType1B = new SType[LDB1 * N1];
1437  OType2B = new SType[LDB2 * N2];
1438 
1439  SIDE = RandomSIDE();
1440  UPLO = RandomUPLO();
1441  TRANSA = RandomTRANS();
1442  DIAG = RandomDIAG();
1443 
1444  if(Teuchos::ESideChar[SIDE] == 'L') // The operator is on the left side
1445  {
1446  LDA1 = GetRandom(MVMIN, MVMAX);
1447  while (LDA1 < M1) {
1448  LDA1 = GetRandom(MVMIN, MVMAX);
1449  }
1450  LDA2 = ConvertType( LDA1, LDA2 );
1451 
1452  OType1A = new SType[LDA1 * M1];
1453  OType2A = new SType[LDA2 * M2];
1454 
1455  for(j1 = 0, j2 = 0; j1 < M1; j1++, j2++)
1456  {
1457  if(Teuchos::EUploChar[UPLO] == 'U') {
1458  // The operator is upper triangular, make sure that the entries are
1459  // only in the upper triangular part of A and the diagonal is non-zero.
1460  for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++)
1461  {
1462  if(k1 < j1) {
1463  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1464  } else {
1465  OType1A[j1*LDA1+k1] = STypezero;
1466  }
1467  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1468  if(k1 == j1) {
1469  if (Teuchos::EDiagChar[DIAG] == 'N') {
1470  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1471  while (OType1A[j1*LDA1+k1] == STypezero) {
1472  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1473  }
1474  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1475  } else {
1476  OType1A[j1*LDA1+k1] = STypeone;
1477  OType2A[j2*LDA2+k2] = STypeone;
1478  }
1479  }
1480  }
1481  } else {
1482  // The operator is lower triangular, make sure that the entries are
1483  // only in the lower triangular part of A and the diagonal is non-zero.
1484  for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++)
1485  {
1486  if(k1 > j1) {
1487  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1488  } else {
1489  OType1A[j1*LDA1+k1] = STypezero;
1490  }
1491  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1492  if(k1 == j1) {
1493  if (Teuchos::EDiagChar[DIAG] == 'N') {
1494  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1495  while (OType1A[j1*LDA1+k1] == STypezero) {
1496  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1497  }
1498  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1499  } else {
1500  OType1A[j1*LDA1+k1] = STypeone;
1501  OType2A[j2*LDA2+k2] = STypeone;
1502  }
1503  }
1504  } // end for(k=0 ...
1505  } // end if(UPLO == 'U') ...
1506  } // end for(j=0 ...
1507  } // if(SIDE == 'L') ...
1508  else // The operator is on the right side
1509  {
1510  LDA1 = GetRandom(MVMIN, MVMAX);
1511  while (LDA1 < N1) {
1512  LDA1 = GetRandom(MVMIN, MVMAX);
1513  }
1514  LDA2 = ConvertType( LDA1, LDA2 );
1515 
1516  OType1A = new SType[LDA1 * N1];
1517  OType2A = new SType[LDA2 * N2];
1518 
1519  for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++)
1520  {
1521  if(Teuchos::EUploChar[UPLO] == 'U') {
1522  // The operator is upper triangular, make sure that the entries are
1523  // only in the upper triangular part of A and the diagonal is non-zero.
1524  for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++)
1525  {
1526  if(k1 < j1) {
1527  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1528  } else {
1529  OType1A[j1*LDA1+k1] = STypezero;
1530  }
1531  OType2A[j1*LDA1+k1] = OType1A[j1*LDA1+k1];
1532  if(k1 == j1) {
1533  if (Teuchos::EDiagChar[DIAG] == 'N') {
1534  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1535  while (OType1A[j1*LDA1+k1] == STypezero) {
1536  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1537  }
1538  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1539  } else {
1540  OType1A[j1*LDA1+k1] = STypeone;
1541  OType2A[j2*LDA2+k2] = STypeone;
1542  }
1543  }
1544  }
1545  } else {
1546  // The operator is lower triangular, make sure that the entries are
1547  // only in the lower triangular part of A and the diagonal is non-zero.
1548  for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++)
1549  {
1550  if(k1 > j1) {
1551  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1552  } else {
1553  OType1A[j1*LDA1+k1] = STypezero;
1554  }
1555  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1556  if(k1 == j1) {
1557  if (Teuchos::EDiagChar[DIAG] == 'N') {
1558  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1559  while (OType1A[j1*LDA1+k1] == STypezero) {
1560  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1561  }
1562  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1563  } else {
1564  OType1A[j1*LDA1+k1] = STypeone;
1565  OType2A[j2*LDA2+k2] = STypeone;
1566  }
1567  }
1568  } // end for(k=0 ...
1569  } // end if(UPLO == 'U') ...
1570  } // end for(j=0 ...
1571  } // end if(SIDE == 'L') ...
1572 
1573  // Fill in the right hand side block B.
1574  for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++) {
1575  for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++) {
1576  OType1B[j1*LDB1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1577  OType2B[j2*LDB2+k2] = OType1B[j1*LDB1+k1];
1578  }
1579  }
1580  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
1581  OType2alpha = OType1alpha;
1582  if(debug)
1583  {
1584  std::cout << "Test #" << TotalTestCount << " -- TRMM -- " << std::endl;
1585  std::cout << "SIDE = " << Teuchos::ESideChar[SIDE] << "\t"
1586  << "UPLO = " << Teuchos::EUploChar[UPLO] << "\t"
1587  << "TRANSA = " << Teuchos::ETranspChar[TRANSA] << "\t"
1588  << "DIAG = " << Teuchos::EDiagChar[DIAG] << std::endl;
1589  std::cout << "OType1alpha = " << OType1alpha << std::endl;
1590  std::cout << "OType2alpha = " << OType2alpha << std::endl;
1591  if(Teuchos::ESideChar[SIDE] == 'L') {
1592  PrintMatrix(OType1A, M1, M1, LDA1, "OType1A", matlab);
1593  PrintMatrix(OType2A, M2, M2, LDA2, "OType2A", matlab);
1594  } else {
1595  PrintMatrix(OType1A, N1, N1, LDA1, "OType1A", matlab);
1596  PrintMatrix(OType2A, N2, N2, LDA2, "OType2A", matlab);
1597  }
1598  PrintMatrix(OType1B, M1, N1, LDB1,"OType1B_before_operation", matlab);
1599  PrintMatrix(OType2B, M2, N2, LDB2,"OType2B_before_operation", matlab);
1600  }
1601  TotalTestCount++;
1602  OType1BLAS.TRMM(SIDE, UPLO, TRANSA, DIAG, M1, N1, OType1alpha, OType1A, LDA1, OType1B, LDB1);
1603  OType2BLAS.TRMM(SIDE, UPLO, TRANSA, DIAG, M2, N2, OType2alpha, OType2A, LDA2, OType2B, LDB2);
1604  if(debug)
1605  {
1606  PrintMatrix(OType1B, M1, N1, LDB1, "OType1B_after_operation", matlab);
1607  PrintMatrix(OType2B, M2, N2, LDB2, "OType2B_after_operation", matlab);
1608  }
1609  if (CompareMatrices(OType1B, M1, N1, LDB1, OType2B, M2, N2, LDB2, TOL)==0)
1610  std::cout << "FAILED TEST!!!!!!" << std::endl;
1611  GoodTestSubcount += CompareMatrices(OType1B, M1, N1, LDB1, OType2B, M2, N2, LDB2, TOL);
1612  delete [] OType1A;
1613  delete [] OType1B;
1614  delete [] OType2A;
1615  delete [] OType2B;
1616  }
1617  GoodTestCount += GoodTestSubcount;
1618  if(verbose || debug) std::cout << "TRMM: " << GoodTestSubcount << " of " << TRMMTESTS << " tests were successful." << std::endl;
1619  if(debug) std::cout << std::endl;
1620  //--------------------------------------------------------------------------------
1621  // End TRMM Tests
1622  //--------------------------------------------------------------------------------
1623 
1624  //--------------------------------------------------------------------------------
1625  // Begin TRSM Tests
1626  //--------------------------------------------------------------------------------
1627  GoodTestSubcount = 0;
1628  for(i = 0; i < TRSMTESTS; i++)
1629  {
1630  M1 = GetRandom(MVMIN, MVMAX);
1631  N1 = GetRandom(MVMIN, MVMAX);
1632  M2 = ConvertType( M1, M2 );
1633  N2 = ConvertType( N1, N2 );
1634 
1635  LDB1 = GetRandom(MVMIN, MVMAX);
1636  while (LDB1 < M1) {
1637  LDB1 = GetRandom(MVMIN, MVMAX);
1638  }
1639  LDB2 = ConvertType( LDB1, LDB2 );
1640 
1641  OType1B = new SType[LDB1 * N1];
1642  OType2B = new SType[LDB2 * N2];
1643 
1644  SIDE = RandomSIDE();
1645  UPLO = RandomUPLO();
1646  TRANSA = RandomTRANS();
1647  // Since the entries are integers, we don't want to use the unit diagonal feature,
1648  // this creates ill-conditioned, nearly-singular matrices.
1649  //DIAG = RandomDIAG();
1650  DIAG = Teuchos::NON_UNIT_DIAG;
1651 
1652  if(Teuchos::ESideChar[SIDE] == 'L') // The operator is on the left side
1653  {
1654  LDA1 = GetRandom(MVMIN, MVMAX);
1655  while (LDA1 < M1) {
1656  LDA1 = GetRandom(MVMIN, MVMAX);
1657  }
1658  LDA2 = ConvertType( LDA1, LDA2 );
1659 
1660  OType1A = new SType[LDA1 * M1];
1661  OType2A = new SType[LDA2 * M2];
1662 
1663  for(j1 = 0, j2 = 0; j1 < M1; j1++, j2++)
1664  {
1665  if(Teuchos::EUploChar[UPLO] == 'U') {
1666  // The operator is upper triangular, make sure that the entries are
1667  // only in the upper triangular part of A and the diagonal is non-zero.
1668  for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++)
1669  {
1670  if(k1 < j1) {
1671  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1672  } else {
1673  OType1A[j1*LDA1+k1] = STypezero;
1674  }
1675  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1676  if(k1 == j1) {
1677  if (Teuchos::EDiagChar[DIAG] == 'N') {
1678  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1679  while (OType1A[j1*LDA1+k1] == STypezero) {
1680  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1681  }
1682  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1683  } else {
1684  OType1A[j1*LDA1+k1] = STypeone;
1685  OType2A[j2*LDA2+k2] = STypeone;
1686  }
1687  }
1688  }
1689  } else {
1690  // The operator is lower triangular, make sure that the entries are
1691  // only in the lower triangular part of A and the diagonal is non-zero.
1692  for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++)
1693  {
1694  if(k1 > j1) {
1695  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1696  } else {
1697  OType1A[j1*LDA1+k1] = STypezero;
1698  }
1699  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1700  if(k1 == j1) {
1701  if (Teuchos::EDiagChar[DIAG] == 'N') {
1702  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1703  while (OType1A[j1*LDA1+k1] == STypezero) {
1704  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1705  }
1706  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1707  } else {
1708  OType1A[j1*LDA1+k1] = STypeone;
1709  OType2A[j2*LDA2+k2] = STypeone;
1710  }
1711  }
1712  } // end for(k=0 ...
1713  } // end if(UPLO == 'U') ...
1714  } // end for(j=0 ...
1715  } // if(SIDE == 'L') ...
1716  else // The operator is on the right side
1717  {
1718  LDA1 = GetRandom(MVMIN, MVMAX);
1719  while (LDA1 < N1) {
1720  LDA1 = GetRandom(MVMIN, MVMAX);
1721  }
1722  LDA2 = ConvertType( LDA1, LDA2 );
1723 
1724  OType1A = new SType[LDA1 * N1];
1725  OType2A = new SType[LDA2 * N2];
1726 
1727  for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++)
1728  {
1729  if(Teuchos::EUploChar[UPLO] == 'U') {
1730  // The operator is upper triangular, make sure that the entries are
1731  // only in the upper triangular part of A and the diagonal is non-zero.
1732  for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++)
1733  {
1734  if(k1 < j1) {
1735  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1736  } else {
1737  OType1A[j1*LDA1+k1] = STypezero;
1738  }
1739  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1740  if(k1 == j1) {
1741  if (Teuchos::EDiagChar[DIAG] == 'N') {
1742  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1743  while (OType1A[j1*LDA1+k1] == STypezero) {
1744  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1745  }
1746  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1747  } else {
1748  OType1A[j1*LDA1+k1] = STypeone;
1749  OType2A[j2*LDA2+k2] = STypeone;
1750  }
1751  }
1752  }
1753  } else {
1754  // The operator is lower triangular, make sure that the entries are
1755  // only in the lower triangular part of A and the diagonal is non-zero.
1756  for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++)
1757  {
1758  if(k1 > j1) {
1759  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1760  } else {
1761  OType1A[j1*LDA1+k1] = STypezero;
1762  }
1763  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1764  if(k1 == j1) {
1765  if (Teuchos::EDiagChar[DIAG] == 'N') {
1766  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1767  while (OType1A[j1*LDA1+k1] == STypezero) {
1768  OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1769  }
1770  OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1];
1771  } else {
1772  OType1A[j1*LDA1+k1] = STypeone;
1773  OType2A[j2*LDA2+k2] = STypeone;
1774  }
1775  }
1776  } // end for(k=0 ...
1777  } // end if(UPLO == 'U') ...
1778  } // end for(j=0 ...
1779  } // end if(SIDE == 'L') ...
1780 
1781  // Fill in the right hand side block B.
1782  for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++)
1783  {
1784  for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++)
1785  {
1786  OType1B[j1*LDB1+k1] = GetRandom(-SCALARMAX, SCALARMAX);
1787  OType2B[j2*LDB2+k2] = OType1B[j1*LDB1+k1];
1788  }
1789  }
1790 
1791  OType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
1792  OType2alpha = OType1alpha;
1793 
1794  if(debug)
1795  {
1796  std::cout << "Test #" << TotalTestCount << " -- TRSM -- " << std::endl;
1797  std::cout << "SIDE = " << Teuchos::ESideChar[SIDE] << "\t"
1798  << "UPLO = " << Teuchos::EUploChar[UPLO] << "\t"
1799  << "TRANSA = " << Teuchos::ETranspChar[TRANSA] << "\t"
1800  << "DIAG = " << Teuchos::EDiagChar[DIAG] << std::endl;
1801  std::cout << "M1="<<M1 << "\t" << "N1="<<N1 << "\t" << "LDA1="<<LDA1 << "\t" << "LDB1="<<LDB1 << std::endl;
1802  std::cout << "M2="<<M2 << "\t" << "N2="<<N2 << "\t" << "LDA2="<<LDA2 << "\t" << "LDB2="<<LDB2 << std::endl;
1803  std::cout << "OType1alpha = " << OType1alpha << std::endl;
1804  std::cout << "OType2alpha = " << OType2alpha << std::endl;
1805  if (Teuchos::ESideChar[SIDE] == 'L') {
1806  PrintMatrix(OType1A, M1, M1, LDA1, "OType1A", matlab);
1807  PrintMatrix(OType2A, M2, M2, LDA2, "OType2A", matlab);
1808  } else {
1809  PrintMatrix(OType1A, N1, N1, LDA1, "OType1A", matlab);
1810  PrintMatrix(OType2A, N2, N2, LDA2, "OType2A", matlab);
1811  }
1812  PrintMatrix(OType1B, M1, N1, LDB1, "OType1B_before_operation", matlab);
1813  PrintMatrix(OType2B, M2, N2, LDB2, "OType2B_before_operation", matlab);
1814  }
1815  TotalTestCount++;
1816 
1817  OType1BLAS.TRSM(SIDE, UPLO, TRANSA, DIAG, M1, N1, OType1alpha, OType1A, LDA1, OType1B, LDB1);
1818  OType2BLAS.TRSM(SIDE, UPLO, TRANSA, DIAG, M2, N2, OType2alpha, OType2A, LDA2, OType2B, LDB2);
1819 
1820  if(debug)
1821  {
1822  PrintMatrix(OType1B, M1, N1, LDB1, "OType1B_after_operation", matlab);
1823  PrintMatrix(OType2B, M2, N2, LDB2, "OType2B_after_operation", matlab);
1824  }
1825 
1826  if (CompareMatrices(OType1B, M1, N1, LDB1, OType2B, M2, N2, LDB2, TOL)==0)
1827  std::cout << "FAILED TEST!!!!!!" << std::endl;
1828  GoodTestSubcount += CompareMatrices(OType1B, M1, N1, LDB1, OType2B, M2, N2, LDB2, TOL);
1829 
1830  delete [] OType1A;
1831  delete [] OType1B;
1832  delete [] OType2A;
1833  delete [] OType2B;
1834  }
1835  GoodTestCount += GoodTestSubcount;
1836  if(verbose || debug) std::cout << "TRSM: " << GoodTestSubcount << " of " << TRSMTESTS << " tests were successful." << std::endl;
1837  if(debug) std::cout << std::endl;
1838  //--------------------------------------------------------------------------------
1839  // End TRSM Tests
1840  //--------------------------------------------------------------------------------
1841 
1842  if((((TotalTestCount - 1) - GoodTestCount) != 0) || (verbose) || (debug))
1843  {
1844  std::cout << GoodTestCount << " of " << (TotalTestCount - 1) << " total tests were successful." << std::endl;
1845  }
1846 
1847  if ((TotalTestCount-1) == GoodTestCount) {
1848  std::cout << "End Result: TEST PASSED" << std::endl;
1849  return 0;
1850  }
1851 
1852  std::cout << "End Result: TEST FAILED" << std::endl;
1853  return (TotalTestCount-GoodTestCount-1);
1854 }
1855 
1856 template<typename TYPE>
1857 TYPE GetRandom(TYPE Low, TYPE High)
1858 {
1859  return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
1860 }
1861 
1862 template<typename T>
1863 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
1864 {
1865  T lowMag = Low.real();
1866  T highMag = High.real();
1867  T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
1868  T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
1869  return std::complex<T>( real, imag );
1870 }
1871 
1872 template<>
1873 int GetRandom(int Low, int High)
1874 {
1875  return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
1876 }
1877 
1878 template<>
1879 double GetRandom(double Low, double High)
1880 {
1881  return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
1882 }
1883 
1884 template<typename TYPE, typename OTYPE>
1885 void PrintVector(TYPE* Vector, OTYPE Size, std::string Name, bool Matlab)
1886 {
1887  std::cout << Name << " =" << std::endl;
1888  OTYPE i;
1889  if(Matlab) std::cout << "[";
1890  for(i = 0; i < Size; i++)
1891  {
1892  std::cout << Vector[i] << " ";
1893  }
1894  if(Matlab) std::cout << "]";
1895  if(!Matlab)
1896  {
1897  std::cout << std::endl << std::endl;
1898  }
1899  else
1900  {
1901  std::cout << ";" << std::endl;
1902  }
1903 }
1904 
1905 template<typename TYPE, typename OTYPE>
1906 void PrintMatrix(TYPE* Matrix, OTYPE Rows, OTYPE Columns, OTYPE LDM, std::string Name, bool Matlab)
1907 {
1908  if(!Matlab)
1909  {
1910  std::cout << Name << " =" << std::endl;
1911  OTYPE i, j;
1912  for(i = 0; i < Rows; i++)
1913  {
1914  for(j = 0; j < Columns; j++)
1915  {
1916  std::cout << Matrix[i + (j * LDM)] << " ";
1917  }
1918  std::cout << std::endl;
1919  }
1920  std::cout << std::endl;
1921  }
1922  else
1923  {
1924  std::cout << Name << " = ";
1925  OTYPE i, j;
1926  std::cout << "[";
1927  for(i = 0; i < Rows; i++)
1928  {
1929  std::cout << "[";
1930  for(j = 0; j < Columns; j++)
1931  {
1932  std::cout << Matrix[i + (j * LDM)] << " ";
1933  }
1934  std::cout << "];";
1935  }
1936  std::cout << "];" << std::endl;
1937  }
1938 }
1939 
1940 template<typename TYPE>
1941 bool CompareScalars(TYPE Scalar1, TYPE Scalar2, typename ScalarTraits<TYPE>::magnitudeType Tolerance)
1942 {
1943  typename ScalarTraits<TYPE>::magnitudeType temp = ScalarTraits<TYPE>::magnitude(Scalar2);
1944  typename ScalarTraits<TYPE>::magnitudeType temp2 = ScalarTraits<TYPE>::magnitude(Scalar1 - Scalar2);
1945  if (temp != ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero()) {
1946  temp2 /= temp;
1947  }
1948  return( temp2 < Tolerance );
1949 }
1950 
1951 
1952 /* Function: CompareVectors
1953  Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
1954 */
1955 template<typename TYPE, typename OTYPE1, typename OTYPE2>
1956 bool CompareVectors(TYPE* Vector1, OTYPE1 Size1, TYPE* Vector2, OTYPE2 Size2, typename ScalarTraits<TYPE>::magnitudeType Tolerance)
1957 {
1958  TYPE temp = ScalarTraits<TYPE>::zero();
1959  typename ScalarTraits<TYPE>::magnitudeType temp2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero();
1960  typename ScalarTraits<TYPE>::magnitudeType temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero();
1961  typename ScalarTraits<TYPE>::magnitudeType sum = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero();
1962  typename ScalarTraits<TYPE>::magnitudeType sum2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero();
1963  OTYPE1 i1;
1964  OTYPE2 i2;
1965  for(i1 = 0, i2 = 0; i1 < Size1; i1++, i2++)
1966  {
1967  sum2 += ScalarTraits<TYPE>::magnitude(ScalarTraits<TYPE>::conjugate(Vector2[i2])*Vector2[i2]);
1968  temp = Vector1[i1] - Vector2[i2];
1969  sum += ScalarTraits<TYPE>::magnitude(ScalarTraits<TYPE>::conjugate(temp)*temp);
1970  }
1971  temp2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum2);
1972  if (temp2 != ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero())
1973  temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum)/temp2;
1974  else
1975  temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum);
1976  if (temp3 > Tolerance )
1977  return false;
1978  else
1979  return true;
1980 }
1981 
1982 /* Function: CompareMatrices
1983  Purpose: Compares the difference between two matrices using relative frobenius-norms, i.e. ||M_1-M_2||_F/||M_2||_F
1984 */
1985 template<typename TYPE, typename OTYPE1, typename OTYPE2>
1986 bool CompareMatrices(TYPE* Matrix1, OTYPE1 Rows1, OTYPE1 Columns1, OTYPE1 LDM1,
1987  TYPE* Matrix2, OTYPE2 Rows2, OTYPE2 Columns2, OTYPE2 LDM2,
1988  typename ScalarTraits<TYPE>::magnitudeType Tolerance)
1989 {
1990  TYPE temp = ScalarTraits<TYPE>::zero();
1991  typename ScalarTraits<TYPE>::magnitudeType temp2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero();
1992  typename ScalarTraits<TYPE>::magnitudeType temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero();
1993  typename ScalarTraits<TYPE>::magnitudeType sum = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero();
1994  typename ScalarTraits<TYPE>::magnitudeType sum2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero();
1995  OTYPE1 i1, j1;
1996  OTYPE2 i2, j2;
1997  for(j1 = 0, j2 = 0; j1 < Columns1; j1++, j2++)
1998  {
1999  for(i1 = 0, i2 = 0; i1 < Rows1; i1++, i2++)
2000  {
2001  sum2 = ScalarTraits<TYPE>::magnitude(ScalarTraits<TYPE>::conjugate(Matrix2[j2*LDM2 + i2])*Matrix2[j2*LDM2 + i2]);
2002  temp = Matrix1[j1*LDM1 + i1] - Matrix2[j2*LDM2 + i2];
2003  sum = ScalarTraits<TYPE>::magnitude(ScalarTraits<TYPE>::conjugate(temp)*temp);
2004  }
2005  }
2006  temp2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum2);
2007  if (temp2 != ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero())
2008  temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum)/temp2;
2009  else
2010  temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum);
2011  if (temp3 > Tolerance)
2012  return false;
2013  else
2014  return true;
2015 }
2016 
2017 
2018 Teuchos::ESide RandomSIDE()
2019 {
2020  Teuchos::ESide result;
2021  int r = GetRandom(1, 2);
2022  if(r == 1)
2023  {
2024  result = Teuchos::LEFT_SIDE;
2025  }
2026  else
2027  {
2028  result = Teuchos::RIGHT_SIDE;
2029  }
2030  return result;
2031 }
2032 
2033 Teuchos::EUplo RandomUPLO()
2034 {
2035  Teuchos::EUplo result;
2036  int r = GetRandom(1, 2);
2037  if(r == 1)
2038  {
2039  result = Teuchos::UPPER_TRI;
2040  }
2041  else
2042  {
2043  result = Teuchos::LOWER_TRI;
2044  }
2045  return result;
2046 }
2047 
2048 Teuchos::ETransp RandomTRANS()
2049 {
2050  Teuchos::ETransp result;
2051  int r = GetRandom(1, 4);
2052  if(r == 1 || r == 2)
2053  {
2054  result = Teuchos::NO_TRANS;
2055  }
2056  else if(r == 3)
2057  {
2058  result = Teuchos::TRANS;
2059  }
2060  else
2061  {
2062  result = Teuchos::CONJ_TRANS;
2063  }
2064  return result;
2065 }
2066 
2067 Teuchos::EDiag RandomDIAG()
2068 {
2069  Teuchos::EDiag result;
2070  int r = GetRandom(1, 2);
2071  if(r == 1)
2072  {
2073  result = Teuchos::NON_UNIT_DIAG;
2074  }
2075  else
2076  {
2077  result = Teuchos::UNIT_DIAG;
2078  }
2079  return result;
2080 }
2081 
SYMMTESTS
#define SYMMTESTS
Definition: cxx_tmpl_main_comp.cpp:71
ROTGTESTS
#define ROTGTESTS
Definition: cxx_tmpl_main_comp.cpp:58
Teuchos::DefaultBLASImpl::TRSM
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const
Solves the matrix equations: op(A)*X=alpha*B or X*op(A)=alpha*B where X and B are m by n matrices...
Definition: Teuchos_BLAS.hpp:1930
Teuchos::DefaultBLASImpl::GER
void GER(const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, const y_type *y, const OrdinalType &incy, ScalarType *A, const OrdinalType &lda) const
Performs the rank 1 operation: A &lt;- alpha*x*y&#39;+A.
Definition: Teuchos_BLAS.hpp:1034
Teuchos::DefaultBLASImpl::AXPY
void AXPY(const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
Perform the operation: y &lt;- y+alpha*x.
Definition: Teuchos_BLAS.hpp:561
Teuchos::DefaultBLASImpl::TRMV
void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType &n, const A_type *A, const OrdinalType &lda, ScalarType *x, const OrdinalType &incx) const
Performs the matrix-vector operation: x &lt;- A*x or x &lt;- A&#39;*x where A is a unit/non-unit n by n upper/low...
Definition: Teuchos_BLAS.hpp:839
MVMAX
#define MVMAX
Definition: cxx_tmpl_main_comp.cpp:47
MVMIN
#define MVMIN
Definition: cxx_tmpl_main_comp.cpp:46
Teuchos_BLAS.hpp
Templated interface class to BLAS routines.
Teuchos::DefaultBLASImpl::GEMV
void GEMV(ETransp trans, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const x_type *x, const OrdinalType &incx, const beta_type beta, ScalarType *y, const OrdinalType &incy) const
Performs the matrix-vector operation: y &lt;- alpha*A*x+beta*y or y &lt;- alpha*A&#39;*x+beta*y where A is a gene...
Definition: Teuchos_BLAS.hpp:686
Teuchos::DefaultBLASImpl::ROT
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
Applies a Givens plane rotation.
Definition: Teuchos_BLAS.hpp:495
TRMMTESTS
#define TRMMTESTS
Definition: cxx_tmpl_main_comp.cpp:73
SType
#define SType
Definition: cxx_tmpl_main_comp.cpp:37
GERTESTS
#define GERTESTS
Definition: cxx_tmpl_main_comp.cpp:68
Teuchos_Time.hpp
Basic wall-clock timer class.
OTYPE
#define OTYPE
Definition: numerics/test/DenseMatrix/cxx_main.cpp:15
IAMAXTESTS
#define IAMAXTESTS
Definition: cxx_tmpl_main_comp.cpp:64
Teuchos::BLAS
Templated BLAS wrapper.
Definition: Teuchos_BLAS.hpp:212
Teuchos::DefaultBLASImpl::NRM2
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Compute the 2-norm of the vector x.
Definition: Teuchos_BLAS.hpp:629
GEMMTESTS
#define GEMMTESTS
Definition: cxx_tmpl_main_comp.cpp:70
Teuchos::DefaultBLASImpl::IAMAX
OrdinalType IAMAX(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Return the index of the element of x with the maximum magnitude.
Definition: Teuchos_BLAS.hpp:650
TRMVTESTS
#define TRMVTESTS
Definition: cxx_tmpl_main_comp.cpp:69
Teuchos::DefaultBLASImpl::GEMM
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
General matrix-matrix multiply.
Definition: Teuchos_BLAS.hpp:1109
PrintMatrix
void PrintMatrix(TYPE *Matrix, OType Rows, OType Columns, OType LDM, std::string Name, bool Matlab=0)
Definition: numerics/test/BLAS/cxx_main.cpp:1768
Teuchos::GlobalMPISession
Initialize, finalize, and query the global MPI session.
Definition: Teuchos_GlobalMPISession.hpp:81
Teuchos::EDiagChar
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char EDiagChar[]
Definition: Teuchos_BLAS.cpp:37
Teuchos::ScalarTraits
This structure defines some basic traits for a scalar field type.
Definition: Teuchos_ScalarTraitsDecl.hpp:58
OType2
#define OType2
Definition: cxx_tmpl_main_comp.cpp:42
Teuchos::DefaultBLASImpl::ASUM
ScalarTraits< ScalarType >::magnitudeType ASUM(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Sum the absolute values of the entries of x.
Definition: Teuchos_BLAS.hpp:582
RandomUPLO
Teuchos::EUplo RandomUPLO()
Definition: numerics/test/BLAS/cxx_main.cpp:1924
DOTTESTS
#define DOTTESTS
Definition: cxx_tmpl_main_comp.cpp:63
Teuchos::DefaultBLASImpl::COPY
void COPY(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
Copy the vector x to the vector y.
Definition: Teuchos_BLAS.hpp:540
GEMVTESTS
#define GEMVTESTS
Definition: cxx_tmpl_main_comp.cpp:67
GetRandom
TYPE GetRandom(TYPE, TYPE)
Definition: numerics/test/BLAS/cxx_main.cpp:1729
ASUMTESTS
#define ASUMTESTS
Definition: cxx_tmpl_main_comp.cpp:60
RandomSIDE
Teuchos::ESide RandomSIDE()
Definition: numerics/test/BLAS/cxx_main.cpp:1909
Teuchos::EUploChar
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char EUploChar[]
Definition: Teuchos_BLAS.cpp:36
RandomDIAG
Teuchos::EDiag RandomDIAG()
Definition: numerics/test/BLAS/cxx_main.cpp:1958
SCALARMAX
#define SCALARMAX
Definition: cxx_tmpl_main_comp.cpp:55
Teuchos::DefaultBLASImpl::DOT
ScalarType DOT(const OrdinalType &n, const x_type *x, const OrdinalType &incx, const y_type *y, const OrdinalType &incy) const
Form the dot product of the vectors x and y.
Definition: Teuchos_BLAS.hpp:606
CompareVectors
bool CompareVectors(TYPE1 *Vector1, TYPE2 *Vector2, OType Size, TYPE2 Tolerance)
Definition: numerics/test/BLAS/cxx_main.cpp:1819
CompareScalars
bool CompareScalars(TYPE1 Scalar1, TYPE2 Scalar2, TYPE2 Tolerance)
Definition: numerics/test/BLAS/cxx_main.cpp:1803
Teuchos::ESideChar
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ESideChar[]
Definition: Teuchos_BLAS.cpp:34
ROTTESTS
#define ROTTESTS
Definition: cxx_tmpl_main_comp.cpp:59
Teuchos::Teuchos_Version
std::string Teuchos_Version()
Definition: Teuchos_Version.hpp:22
PrintVector
void PrintVector(TYPE *Vector, OType Size, std::string Name, bool Matlab=0)
Definition: numerics/test/BLAS/cxx_main.cpp:1747
main
int main(int argc, char *argv[])
Definition: core/example/CommandLineProcessor/cxx_main.cpp:21
CompareMatrices
bool CompareMatrices(TYPE1 *Matrix1, TYPE2 *Matrix2, OType Rows, OType Columns, OType LDM, TYPE2 Tolerance)
Definition: numerics/test/BLAS/cxx_main.cpp:1849
ConvertType
TYPE2 ConvertType(TYPE1, TYPE2)
Definition: numerics/test/BLAS/cxx_main.cpp:1880
Teuchos::DefaultBLASImpl::SYMM
void SYMM(ESide side, EUplo uplo, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
Performs the matrix-matrix operation: C &lt;- alpha*A*B+beta*C or C &lt;- alpha*B*A+beta*C where A is an m ...
Definition: Teuchos_BLAS.hpp:1393
Teuchos::DefaultBLASImpl::TRMM
void TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const
Performs the matrix-matrix operation: B &lt;- alpha*op(A)*B or B &lt;- alpha*B*op(A) where op(A) is an unit...
Definition: Teuchos_BLAS.hpp:1697
Teuchos_GlobalMPISession.hpp
A MPI utilities class, providing methods for initializing, finalizing, and querying the global MPI se...
Teuchos::DefaultBLASImpl::SCAL
void SCAL(const OrdinalType &n, const ScalarType &alpha, ScalarType *x, const OrdinalType &incx) const
Scale the vector x by the constant alpha.
Definition: Teuchos_BLAS.hpp:522
Teuchos::DefaultBLASImpl::ROTG
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
Computes a Givens plane rotation.
Definition: Teuchos_BLAS.hpp:485
COPYTESTS
#define COPYTESTS
Definition: cxx_tmpl_main_comp.cpp:62
TRSMTESTS
#define TRSMTESTS
Definition: cxx_tmpl_main_comp.cpp:74
Teuchos::ETranspChar
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
Definition: Teuchos_BLAS.cpp:35
Teuchos::DefaultBLASImpl::SYRK
void SYRK(EUplo uplo, ETransp trans, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
Performs the symmetric rank k operation: C &lt;- alpha*A*A&#39;+beta*C or C &lt;- alpha*A&#39;*A+beta*C, where alpha and beta are scalars, C is an n by n symmetric matrix and A is an n by k matrix in the first case or k by n matrix in the second case.
Definition: Teuchos_BLAS.hpp:1523
RandomTRANS
Teuchos::ETransp RandomTRANS()
Definition: numerics/test/BLAS/cxx_main.cpp:1939
NRM2TESTS
#define NRM2TESTS
Definition: cxx_tmpl_main_comp.cpp:65
AXPYTESTS
#define AXPYTESTS
Definition: cxx_tmpl_main_comp.cpp:61
OType1
#define OType1
Definition: cxx_tmpl_main_comp.cpp:41
SCALTESTS
#define SCALTESTS
Definition: cxx_tmpl_main_comp.cpp:66
SYRKTESTS
#define SYRKTESTS
Definition: cxx_tmpl_main_comp.cpp:72
Teuchos_Version.hpp
Generated on Fri Dec 20 2024 09:21:21 for Teuchos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5