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