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