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