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