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