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