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