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