MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DenseLinAlgPack_TestGenMatrixOp.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) 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 Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 //
42 // 2/10/00: g++ 2.95.2 is giving me some trouble with update_success(...) in template
43 // functions for some reason?
44 //
45 
46 #include <iomanip>
47 #include <ostream>
48 #include <vector>
49 #include <typeinfo>
50 
51 #include "DenseLinAlgPack_TestDenseLinAlgPack.hpp"
52 #include "DenseLinAlgPack_DMatrixClass.hpp"
53 #include "DenseLinAlgPack_DMatrixOp.hpp"
54 #include "DenseLinAlgPack_DVectorOut.hpp"
55 #include "DenseLinAlgPack_DMatrixOut.hpp"
56 #include "DenseLinAlgPack_MatVecCompare.hpp"
57 
58 namespace {
59 
60 using DenseLinAlgPack::size_type;
61 using DenseLinAlgPack::value_type;
62 using DenseLinAlgPack::DVector;
63 using DenseLinAlgPack::DVectorSlice;
64 using DenseLinAlgPack::DMatrix;
65 using DenseLinAlgPack::DMatrixSlice;
66 
67 // vs_lhs = alpha * op(M_rhs1) * vs_rhs2 + beta * vs_lhs
68 template<class M_t>
69 void test_MtV( M_t& M_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2
70  , const DVectorSlice& expected_MtV, DVector* tmp1, DVector* tmp2
71  , std::ostream* out, bool* success )
72 {
73  using BLAS_Cpp::rows;
74  using DenseLinAlgPack::Vt_S;
75  using DenseLinAlgPack::Vp_StV;
76  using DenseLinAlgPack::Vp_StMtV;
77  using DenseLinAlgPack::comp;
78  using TestingHelperPack::update_success;
79 
80  // Check alpha = 1.0, 0.5. Check beta = 0.0, 0.5, 1.0
81  bool result = true;
82  value_type
83  alpha = 1.0,
84  beta = 0.0;
85 
86  const size_type
87  m = M_rhs1.rows(),
88  n = M_rhs1.cols();
89 
90  tmp1->resize( rows(m,n,trans_rhs1) );
91  tmp2->resize( tmp1->dim() );
92 
93  for( int i = 0; i < 2; ++i, alpha -= 0.5 ) {
94  beta = 0.0;
95  for( int j = 0; j < 3; ++j, beta += 0.5 ) {
96  if(out)
97  *out << "vs_lhs = " << alpha << " * M_rhs1"
98  << ( trans_rhs1 == BLAS_Cpp::trans ? '\'' : ' ' )
99  << " * vs_rhs2 + " << beta << " * vs_lhs : ";
100  // Initialize vs_lhs
101  for( int k = 1; k <= tmp1->dim(); ++k )
102  (*tmp1)(k) = (*tmp2)(k) = k;
103  // Compute expected
104  Vt_S( &(*tmp1)(), beta );
105  Vp_StV( &(*tmp1)(), alpha, expected_MtV );
106  // Computed
107  Vp_StMtV( &(*tmp2)(), alpha, M_rhs1, trans_rhs1, vs_rhs2, beta );
108  // Compare
109  update_success( result = comp( *tmp1, *tmp2 ), success );
110  if(out) *out << result << std::endl;
111  if(out && !result) *out << "expected_vs_lhs =\n" << *tmp1
112  << "computed_vs_lhs =\n"<< *tmp2;
113  }
114  }
115 }
116 
117 // Print a '\'' for trans and ' ' for no_trans
118 inline char trans_char(BLAS_Cpp::Transp _trans) {
119  return (_trans == BLAS_Cpp::trans ? '\'' : ' ' );
120 }
121 
122 // Test DenseLinAlgPack compatable Level-3 BLAS (Mp_StMtM(...))
123 template <class M_t>
124 void test_MtM( M_t& B, BLAS_Cpp::Transp trans_B, const DMatrixSlice& I
125  , BLAS_Cpp::Transp trans_I, std::ostream* out
126  , DMatrix* C1, DMatrix* C2, bool* success )
127 {
128  using BLAS_Cpp::no_trans;
129  using BLAS_Cpp::rows;
130  using BLAS_Cpp::cols;
131  using DenseLinAlgPack::assign;
132  using DenseLinAlgPack::Mp_StM;
133  using DenseLinAlgPack::Mp_StMtM;
134  using DenseLinAlgPack::comp;
135  using TestingHelperPack::update_success;
136 
137  value_type alpha = 1.0, beta = 0.0;
138  bool result = true;
139 
140  for( int i = 1; i <= 3; ++i ) {
141  // Set alpha + (beta)(0.5) = 1.0
142  switch(i) {
143  case 1:
144  alpha = 1.0;
145  beta = 0.0;
146  break;
147  case 2:
148  alpha = 0.75;
149  beta = 0.5;
150  break;
151  case 3:
152  alpha = 0.5;
153  beta = 1.0;
154  break;
155  }
156  // Perform tests
157  const char
158  B_trans_chr = (trans_B==no_trans?' ':'\''),
159  I_trans_chr = (trans_I==no_trans?' ':'\'');
160  // expected result
161  const size_type
162  rows_B = rows( B.rows(), B.cols(), trans_B ),
163  cols_B = cols( B.rows(), B.cols(), trans_B );
164  C1->resize(rows_B,cols_B);
165  (*C1) = 0.0;
166  Mp_StM( &(*C1)(), 1.0, B, trans_B );
167  // (left)
168  if(out) *out << "C = " << alpha << " * B" << B_trans_chr
169  << "* I" << I_trans_chr << " + " << beta << " * (0.5) * B" << B_trans_chr
170  << " : ";
171  C2->resize(rows_B,cols_B);
172  (*C2) = 0.0;
173  Mp_StM( &(*C2)(), 0.5, B, trans_B );
174  size_type n = cols( B.rows(), B.cols(), trans_B );
175  Mp_StMtM( &(*C2)(), alpha, B, trans_B, I(1,n,1,n), trans_I, beta );
176  update_success( result = comp( (*C1)(), (*C2)() ), success );
177  if(out) *out << result << std::endl;
178  if( out && !result ) *out << "C1 =\n" << (*C1) << "C2 =\n" << (*C2);
179  // (right)
180  if(out) *out << "C = " << alpha << " * I" << I_trans_chr
181  << "* B" << B_trans_chr << " + " << beta << " * (0.5) * B" << B_trans_chr
182  << " : ";
183  (*C2) = 0.0;
184  Mp_StM( &(*C2)(), 0.5, B, trans_B );
185  n = rows( B.rows(), B.cols(), trans_B );
186  Mp_StMtM( &(*C2)(), alpha, I(1,n,1,n), trans_I, B, trans_B, beta );
187  update_success( result = comp( (*C1)(), (*C2)() ), success );
188  if(out) *out << result << std::endl;
189  if( out && !result ) *out << "C1 =\n" << (*C1) << "C2 =\n" << (*C2);
190 
191  }
192 }
193 
194 } // end namespace
195 
196 bool DenseLinAlgPack::TestingPack::TestGenMatrixOp(std::ostream* out)
197 {
198 
199  using BLAS_Cpp::trans;
200  using BLAS_Cpp::no_trans;
201  using BLAS_Cpp::upper;
202  using BLAS_Cpp::lower;
203  using BLAS_Cpp::unit;
204  using BLAS_Cpp::nonunit;
205  using BLAS_Cpp::trans_to_bool;
206  using DenseLinAlgPack::comp;
207  using DenseLinAlgPack::sqrt_eps;
208 
209  bool success = true;
210  bool result, result1, result2;
211 
212  if(out)
213  *out << "\n**************************************"
214  << "\n*** Testing GenMatrixOp Operations ***"
215  << "\n**************************************\n"
216  << std::boolalpha << std::setw(10);
217 
218  // Arrays for looping through options.
219 
220  BLAS_Cpp::Uplo a_uplo[2] = { lower, upper };
221  char str_uplo[2][10] = { "lower", "upper" };
222 
223  BLAS_Cpp::Diag a_diag[2] = { unit, nonunit };
224  char str_diag[2][10] = { "unit", "nonunit" };
225 
226  BLAS_Cpp::Transp a_trans[2] = { trans, no_trans };
227 // char str_trans[2][10] = { "trans", "no_trans" };
228 
229  try {
230 
231  const size_type
232  m = 6,
233  n = 8;
234 
235  const value_type
236  ptr[m*n] =
237  { 1.1, 2.1, 3.1, 4.1, 5.1, 6.1,
238  1.2, 2.2, 3.2, 4.2, 5.2, 6.2,
239  1.3, 2.3, 3.3, 4.3, 5.3, 6.3,
240  1.4, 2.4, 3.4, 4.4, 5.4, 6.4,
241  1.5, 2.5, 3.5, 4.5, 5.5, 6.5,
242  1.6, 2.6, 3.6, 4.6, 5.6, 6.6,
243  1.7, 2.7, 3.7, 4.7, 5.7, 6.7,
244  1.8, 2.8, 3.8, 4.8, 5.8, 6.8 };
245 
246  DMatrix gm_lhs(m,n);
247  const DMatrixSlice gms_rhs(const_cast<value_type*>(ptr),m*n,m,m,n);
248 
249  // ////////////////////////////////////////////////////////////////////////////////
250  // Test Element-wise Assignment functions not already tested in TestGenMatrixClass
251 
252  if(out) *out << "\n***\n*** Test Element-wise Assignment functions not already "
253  "tested in TestGenMatrixClass(...)\n***\n";
254 
255  // gms_lhs = op(gms_rhs)
256  if(out) *out << "\ngms_lhs = op(gms_rhs)\n"
257  << "\ngm_lhs.resize(n,m); assign( &gm_lhs(), gms_rhs, trans );\n";
258  gm_lhs.resize(n,m);
259  assign( &gm_lhs(), gms_rhs, trans );
260  update_success( result = comp( gm_lhs(), no_trans, gms_rhs, trans ), &success );
261  if(out && !result) *out << "gm_lhs =\n" << gm_lhs();
262  if(out) *out << "gm_lhs == gms_rhs' : " << result << std::endl;
263 
264  // gms_lhs = op(gms_rhs)
265  if(out) *out << "\ngm_lhs = op(gms_rhs)\n"
266  << "\ngm_lhs.resize(1,1); assign( &gm_lhs, gms_rhs, trans );\n";
267  gm_lhs.resize(1,1);
268  assign( &gm_lhs, gms_rhs, trans );
269  update_success( result = comp( gm_lhs(), no_trans, gms_rhs, trans ), &success );
270  if(out && !result) *out << "gm_lhs =\n" << gm_lhs();
271  if(out) *out << "gm_lhs == gms_rhs' : " << result << std::endl;
272 
273  // tri_ele_gms_lhs = alpha (elementwise)
274 
275  if(out) *out << "\ntri_ele_gms_lhs = alpha (elementwise)\n";
276 
277  if(out) *out << "\ngms_lhs.resize(m,m); gm_lhs = 1.0; assign(&nonconst_tri_ele(gm_lhs(),lower),2.0);\n";
278  gm_lhs.resize(m,m);
279  gm_lhs = 1.0;
280  assign(&nonconst_tri_ele(gm_lhs(),lower),2.0);
281  update_success( result1 = comp(tri_ele(gm_lhs(),lower),2.0), &success );
282  update_success( result2 = comp(tri_ele(gm_lhs(1,m-1,2,m),upper),1.0), &success );
283  if(out && (!result1 || !result2) )
284  *out << "gm_lhs =\n" << gm_lhs();
285  if(out) *out << "tri_ele(gm_lhs(),lower) == 2.0 : " << result1 << std::endl
286  << "tri_ele(gm_lhs(1,m-1,2,m),upper) == 1.0 : " << result2 << std::endl;
287 
288  if(out) *out << "\nassign(&tri_ele(gm_lhs(1,m-1,2,m),upper),3.0);\n";
289  assign(&nonconst_tri_ele(gm_lhs(1,m-1,2,m),upper),3.0);
290  update_success( result1 = comp(tri_ele(gm_lhs(),lower),2.0), &success );
291  update_success( result2 = comp(tri_ele(gm_lhs(1,m-1,2,m),upper),3.0), &success );
292  if(out && (!result1 || !result2) )
293  *out << "gm_lhs =\n" << gm_lhs();
294  if(out) *out << "tri_ele(gm_lhs(),lower) == 2.0 : " << result1 << std::endl
295  << "tri_ele(gm_lhs(1,m-1,2,m),upper) == 3.0 : " << result2 << std::endl;
296 
297  // tri_ele_gms_lhs = tri_ele_gms_rhs
298 
299  if(out) *out << "\ntri_ele_gms_lhs = tri_ele_gms_rhs\n"
300  << "\nassign(&tri_ele(gm_lhs(2,m,1,m-1),lower),tri_ele(gms_rhs(1,m-1,2,m),upper));\n";
301  assign(&nonconst_tri_ele(gm_lhs(2,m,1,m-1),lower),tri_ele(gms_rhs(1,m-1,2,m),upper));
302  if(out) *out << "assign(&tri_ele(gm_lhs(1,m-1,2,m),upper),tri_ele(gms_rhs(2,m,1,m-1),lower));\n";
303  assign(&nonconst_tri_ele(gm_lhs(1,m-1,2,m),upper),tri_ele(gms_rhs(2,m,1,m-1),lower));
304  if(out) *out << "gm_lhs.diag() = gms_rhs.diag();\n";
305  gm_lhs.diag() = gms_rhs.diag();
306  update_success( result1 = comp( gm_lhs(1,m,1,m), no_trans, gms_rhs(1,m,1,m), trans ), &success );
307  update_success( result2 = comp( tri_ele(gm_lhs(2,m,1,m-1),lower), tri_ele(gms_rhs(1,m-1,2,m),upper) ), &success );
308  if(out && (!result1 || !result2) )
309  *out << "gm_lhs =\n" << gm_lhs();
310  if(out) *out << "gm_lhs(1,m,1,m) == gms_rhs(1,m,1,m)' : " << result1 << std::endl
311  << "tri_ele(gm_lhs(2,m,1,m-1),lower) == tri_ele(gms_rhs(1,m-1,2,m),upper) : " << result2 << std::endl;
312 
313  // //////////////////////////////////////////
314  // Test Element-wise Algebraic Operations
315 
316  if(out) *out << "\n***\n*** Test Element-wise Algebraic Operations\n***\n";
317 
318  // gms_lhs *= alpha
319  if(out) *out << "\ngms_lhs *= alpha\n"
320  << "\ngm_lhs = 1.0; Mt_S( &gm_lhs(), 2.0 );\n";
321  gm_lhs = 1.0;
322  Mt_S( &gm_lhs(), 2.0 );
323  update_success( result = comp( gm_lhs(), 2.0 ), &success );
324  if(out && !result)
325  *out << "gm_lhs =\n" << gm_lhs();
326  if(out) *out << "gm_lhs == 2.0: " << result << std::endl;
327 
328  // tri_lhs *= alpha
329  {
330  if(out) *out << "\ntri_lhs *= alpha\n"
331  << "\ngm_lhs = 1.0;\nLet tm1 = tri_ele(gm_lhs(1,m,1,m),BLAS_Cpp::lower), "
332  "tm2 = tri_ele(gm_lhs(1,m-1,2,m),BLAS_Cpp::upper)\n";
333  gm_lhs = 1.0;
334  DMatrixSliceTriEle
335  tm1 = nonconst_tri_ele(gm_lhs(1,m,1,m),BLAS_Cpp::lower),
336  tm2 = nonconst_tri_ele(gm_lhs(1,m-1,2,m),BLAS_Cpp::upper);
337  if(out) *out << "Mt_S( &tm1, 2.0 );\n";
338  Mt_S( &tm1, 2.0 );
339  if(out) *out << "Mt_S( &tm2, 3.0 );\n";
340  Mt_S( &tm2, 3.0 );
341  update_success( result1 = comp( tm1, 2.0 ), &success );
342  update_success( result2 = comp( tm2, 3.0 ), &success );
343  if(out && (!result1 || !result2) )
344  *out << "gm_lhs =\n" << gm_lhs();
345  if(out) *out << "tm1 == 2.0 : " << result1 << std::endl
346  << "tm2 == 3.0 : " << result2 << std::endl;
347  }
348 
349  // tri_lhs += alpha * tri_rhs
350  if(out) *out << "\ntri_lhs += alpha * tri_rhs\n"
351  "gm_lhs = 1.0;\n";
352  gm_lhs = 1.0;
353  if(out) *out << "Mp_StM( &tri_ele(gm_lhs(2,m,1,m-1),lower), 2.0, tri_ele(gm_lhs(1,m-1,2,m),upper) );\n";
354  Mp_StM( &nonconst_tri_ele(gm_lhs(2,m,1,m-1),lower), 2.0, tri_ele(gm_lhs(1,m-1,2,m),upper) );
355  update_success( result1 = comp( tri_ele(gm_lhs(2,m,1,m-1),lower), 3.0 ), &success );
356  update_success( result2 = comp( tri_ele(gm_lhs(1,m,1,m),upper), 1.0 ), &success );
357  if(out && (!result1 || !result2) )
358  *out << "gm_lhs =\n" << gm_lhs();
359  if(out) *out << "tri_ele(gm_lhs(2,m,1,m-1),lower) == 3.0 : " << result1 << std::endl
360  << "tri_ele(gm_lhs(1,m,1,m),upper) == 1.0 : " << result2 << std::endl;
361 
362  // LinAlgOpPack compatable (Mp_StM(...)).
363  if(out) *out << "\n****** LinAlgOpPack compatable (Mp_StM(...))\n";
364 
365  // gms_lhs += alpha * op(gms_rhs)
366  if(out) *out << "\ngms_lhs += alpha * op(gms_rhs)\n"
367  "--------------------------\n"
368  "gm_lhs.resize(m,n); gm_lhs = 0.0;"
369  "Mp_StM( &gm_lhs(), 0.5, gm_rhs(), no_trans ); "
370  "Mp_StM( &gm_lhs(), 0.5, gm_rhs(), no_trans );\n";
371  gm_lhs.resize(m,n);
372  gm_lhs = 0.0;
373  Mp_StM( &gm_lhs(), 0.5, gms_rhs, no_trans );
374  Mp_StM( &gm_lhs(), 0.5, gms_rhs, no_trans );
375  update_success( result = comp( gm_lhs(), gms_rhs ), &success );
376  if(out && !result)
377  *out << "gm_lhs =\n" << gm_lhs();
378  if(out) *out << "gm_lhs == gm_rhs : " << result << std::endl;
379 
380  // gms_lhs += alpha * op(sym_rhs)
381  if(out) *out << "\ngms_lhs += alpha * op(sym_rhs)\n"
382  "------------------------------\n";
383  gm_lhs.resize(m,m);
384  {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
385  for(int i_trans = 0; i_trans < 2; ++i_trans) {
386  const BLAS_Cpp::Uplo _uplo = a_uplo[i_uplo];
387  const BLAS_Cpp::Transp _trans = a_trans[i_trans];
388  if(out)
389  *out
390  << "gms_lhs += alpha * sym(gms_rhs(1,m,1,m),"
391  << str_uplo[i_uplo] << ")" << (_trans == trans ? '\'' : ' ' ) << " : ";
392  gm_lhs = 0.0;
393  const DMatrixSliceSym M = sym(gms_rhs(1,m,1,m),_uplo);
394  Mp_StM( &gm_lhs(), 0.5, M, _trans );
395  Mp_StM( &gm_lhs(), 0.5, M, _trans );
396  update_success( result1 = comp( tri_ele(gm_lhs(),lower), tri_ele(M.gms(),_uplo) )
397  , &success );
398  update_success( result2 = comp( tri_ele(gm_lhs(),upper), tri_ele(M.gms(),_uplo) )
399  , &success );
400  if(out) *out << ( result1 && result2 ) << std::endl;
401  if( out && ( !result1 || !result2 ) )
402  *out << "gm_lhs =\n" << gm_lhs();
403  }
404  }}
405 
406  // gms_lhs += alpha * op(tri_rhs)
407  if(out) *out << "\ngms_lhs += alpha * op(tri_rhs)\n"
408  "------------------------------\n";
409  gm_lhs.resize(m,m);
410  {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
411  for(int i_diag = 0; i_diag < 2; ++i_diag) {
412  for(int i_trans = 0; i_trans < 2; ++i_trans) {
413  const BLAS_Cpp::Uplo _uplo = a_uplo[i_uplo];
414  const BLAS_Cpp::Diag _diag = a_diag[i_diag];
415  const BLAS_Cpp::Transp _trans = a_trans[i_trans];
416  if(out)
417  *out
418  << "gms_lhs += alpha * tri(gms_rhs(1,m,1,m),"
419  << str_uplo[i_uplo] << "," << str_diag[i_diag] << ")"
420  << (_trans == trans ? '\'' : ' ' ) << " : ";
421  // compute
422  gm_lhs = 0.0;
423  const DMatrixSliceTri M = tri(gms_rhs(1,m,1,m),_uplo,_diag);
424  Mp_StM( &gm_lhs(), 0.5, M, _trans );
425  Mp_StM( &gm_lhs(), 0.5, M, _trans );
426  // test diagonal
427  if( _diag == nonunit )
428  result = comp( gm_lhs.diag(), gms_rhs.diag() );
429  else
430  result = comp( gm_lhs.diag(), 1.0 );
431  update_success( result, &success );
432  // test the rest
433  const BLAS_Cpp::Uplo
434  as_uplo = ( ( _uplo == lower && _trans == no_trans )
435  || ( _uplo == upper && _trans == trans )
436  ? lower : upper ),
437  oth_as_uplo = ( as_uplo == lower ? upper : lower );
438  const DMatrixSliceTriEle
439  M_ele = tri_ele( ( _uplo==lower ? M.gms()(2,m,1,m-1) : M.gms()(1,m-1,2,m) )
440  , _uplo ),
441  tri_reg_ele = tri_ele( ( as_uplo==lower ? gm_lhs(2,m,1,m-1) : gm_lhs(1,m-1,2,m) )
442  , as_uplo ),
443  oth_tri_reg_ele = tri_ele( ( oth_as_uplo==lower ? gm_lhs(2,m,1,m-1) : gm_lhs(1,m-1,2,m) )
444  , oth_as_uplo );
445  update_success( result1 = comp( tri_reg_ele, M_ele ), &success );
446  update_success( result2 = comp( oth_tri_reg_ele, 0.0 ), &success );
447  if(out) *out << ( result && result1 && result2 ) << std::endl;
448  if( out && ( !result || !result1 || !result2 ) )
449  *out << "gm_lhs =\n" << gm_lhs();
450  }
451  }
452  }}
453 
454  // //////////////////////////////////////////
455  // Level-2 BLAS
456 
457  if(out) *out << "\n***\n*** Level-2 BLAS \n***\n";
458 
459  DVector tmp1, tmp2, vs_rhs(n), expected_MtV;
460  {for(int k = 1; k <= n; ++k) vs_rhs(k) = k; }
461 
462  // *** Triagnular matrices
463  if(out) *out << "\n*** BLAS-2 Equivalent Triangular Matrix operations\n";
464 
465  {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
466  for(int i_diag = 0; i_diag < 2; ++i_diag) {
467  for(int i_trans = 0; i_trans < 2; ++i_trans) {
468  if(out)
469  *out << "\nLet M = tri(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
470  << "," << str_diag[i_diag] << ")"
471  << (a_trans[i_trans] == trans ? '\'' : ' ' ) << std::endl;
472  DMatrixSliceTri M = tri(gms_rhs(1,m,1,m),a_uplo[i_uplo],a_diag[i_diag]);
473  // v_lhs
474  tmp1.resize(1);
475  tmp2.resize(1);
476  V_InvMtV( &tmp1, M, a_trans[i_trans], vs_rhs(1,m) );
477  V_MtV( &tmp2, M, a_trans[i_trans], tmp1() );
478  update_success( result = comp( tmp2, vs_rhs(1,m)), &success );
479  if(out) *out << "(v_lhs)... M * (inv(M)*vs_rhs(1,m)) == vs_rhs(1,m) : "
480  << result << std::endl;
481  // vs_lhs
482  V_InvMtV( &tmp1(), M, a_trans[i_trans], vs_rhs(1,m) );
483  V_MtV( &tmp1(), M, a_trans[i_trans], tmp1() );
484  update_success( result = comp( tmp1, vs_rhs(1,m)), &success );
485  if(out) *out << "(vs_lhs)... M * (inv(M)*vs_rhs(1,m)) == vs_rhs(1,m) : "
486  << result << std::endl;
487  }
488  }
489  }}
490 
491  // *** LinAlgOpPack complient
492  if(out) *out << "\n*** Test DenseLinAlgPack complient (Vp_StMtV(...))\n";
493 
494  // vs_lhs = alpha * op(gms_rhs1) * vs_rhs2 + beta * vs_lhs (xGEMV)
495 
496  if(out) *out << "\n****** General Matrix MtV\n";
497 
498  // no_trans
499  if(out) *out << "\nvs_lhs = alpha * gms_rhs1 * vs_rhs2 + beta * vs_lhs (xGEMV)\n";
500  expected_MtV.resize(m);
501  expected_MtV = 0.0;
502  {for(int i=1;i<=m;++i)
503  for(int j=1;j<=n;++j)
504  expected_MtV(i) += gms_rhs(i,j) * vs_rhs(j);
505  }
506  test_MtV( gms_rhs, no_trans, vs_rhs(1,n), expected_MtV, &tmp1, &tmp2
507  , out, &success );
508 
509  // trans
510  if(out) *out << "\nvs_lhs = alpha * gms_rhs1' * vs_rhs2 + beta * vs_lhs (xGEMV)\n";
511  expected_MtV.resize(n);
512  expected_MtV = 0.0;
513  {for(int i=1;i<=n;++i)
514  for(int j=1;j<=m;++j)
515  expected_MtV(i) += gms_rhs(j,i) * vs_rhs(j);
516  }
517  test_MtV( gms_rhs, trans, vs_rhs(1,m), expected_MtV, &tmp1, &tmp2
518  , out, &success );
519 
520  // vs_lhs = alpha * op(sym_rhs1) * vs_rhs2 + beta * vs_lhs (xSYMV)
521 
522  if(out) *out << "\n****** Symmetric Matrix MtV\n";
523 
524  {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
525  for(int i_trans = 0; i_trans < 2; ++i_trans) {
526  if(out)
527  *out
528  << "\nvs_lhs = alpha * sym(gms_rhs(1,m,1,m),"
529  << str_uplo[i_uplo] << ")"
530  << (a_trans[i_trans] == trans ? '\'' : ' ' )
531  << " * vs_rhs2 + beta * vs_lhs (xSYMV)\n";
532  // compute expected MtV
533  expected_MtV.resize(m);
534  expected_MtV = 0.0;
535  switch(a_uplo[i_uplo]) {
536  case lower: {
537  for(int i=1;i<=m;++i)
538  for(int j=1;j<=m;++j)
539  expected_MtV(i)
540  += (j < i ? gms_rhs(i,j) : gms_rhs(j,i) ) * vs_rhs(j);
541  break;
542  }
543  case upper: {
544  for(int i=1;i<=m;++i)
545  for(int j=1;j<=m;++j)
546  expected_MtV(i)
547  += (j > i ? gms_rhs(i,j) : gms_rhs(j,i) ) * vs_rhs(j);
548  break;
549  }
550  }
551  test_MtV( sym(gms_rhs(1,m,1,m),a_uplo[i_uplo]), a_trans[i_trans]
552  , vs_rhs(1,m), expected_MtV, &tmp1, &tmp2, out, &success );
553  }
554  }}
555 
556 
557  // vs_lhs = alpha * op(tri_rhs1) * vs_rhs2 + beta * vs_lhs.
558 
559  if(out) *out << "\n****** Triangular Matrix MtV\n";
560 
561  {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
562  for(int i_diag = 0; i_diag < 2; ++i_diag) {
563  for(int i_trans = 0; i_trans < 2; ++i_trans) {
564  const BLAS_Cpp::Uplo _uplo = a_uplo[i_uplo];
565  const BLAS_Cpp::Diag _diag = a_diag[i_diag];
566  const BLAS_Cpp::Transp _trans = a_trans[i_trans];
567  if(out)
568  *out
569  << "\nvs_lhs = alpha * tri(gms_rhs(1,m,1,m),"
570  << str_uplo[i_uplo] << "," << str_diag[i_diag]
571  << ")" << (_trans == trans ? '\'' : ' ' )
572  << " * vs_rhs2 + beta * vs_lhs\n";
573  DMatrixSliceTri M = tri(gms_rhs(1,m,1,m),_uplo,_diag);
574  //
575  // Compute expected MtV
576  //
577  // Determine conceptually lower or upper triangular
578  const BLAS_Cpp::Uplo
579  as_uplo = ( ( _uplo == lower && _trans == no_trans )
580  || ( _uplo == upper && _trans == trans )
581  ? lower : upper );
582  // Compute expected
583  expected_MtV.resize(m);
584  expected_MtV = 0.0;
585  switch(as_uplo) {
586  case lower: {
587  for(int i=1;i<=m;++i) {
588  for(int j=1;j<i;++j) {
589  const int
590  _i = ( _trans == no_trans ? i : j ),
591  _j = ( _trans == no_trans ? j : i );
592  expected_MtV(i) += gms_rhs(_i,_j) * vs_rhs(j);
593  }
594  expected_MtV(i) += (_diag==unit?1.0:gms_rhs(i,i))*vs_rhs(i);
595  }
596  break;
597  }
598  case upper: {
599  for(int i=1;i<=m;++i) {
600  expected_MtV(i) += (_diag==unit?1.0:gms_rhs(i,i))*vs_rhs(i);
601  for(int j=i+1;j<=m;++j) {
602  const int
603  _i = ( _trans == no_trans ? i : j ),
604  _j = ( _trans == no_trans ? j : i );
605  expected_MtV(i) += gms_rhs(_i,_j) * vs_rhs(j);
606  }
607  }
608  break;
609  }
610  }
611  test_MtV( tri(gms_rhs(1,m,1,m),_uplo,_diag), _trans, vs_rhs(1,m)
612  , expected_MtV, &tmp1, &tmp2, out, &success );
613  }
614  }
615  }}
616 
617  // *** Symmetric Matrix Updates.
618 
619  if(out) *out << "\n*** Symmetric Matrix Updates\n";
620 
621  // sym_lhs = alpha * vs_rhs * vs_rhs' + sym_lhs (BLAS xSYR).
622  if(out) *out << "\n****** sym_lhs = alpha * vs_rhs * vs_rhs' + sym_lhs (BLAS xSYR)\n";
623 
624  gm_lhs.resize(m,m);
625  gm_lhs = 0.0;
626  DMatrix ex_gm_lhs(0.0,m,m);
627  value_type alpha = 2.0;
628 
629  {for( int i_uplo = 0; i_uplo < 2; ++i_uplo) {
630  const BLAS_Cpp::Uplo _uplo = a_uplo[i_uplo];
631  if(out)
632  *out << "\nFor sym_lhs = sym(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
633  << ")\n";
634  assign( &nonconst_tri_ele(gm_lhs(),_uplo), tri_ele(gms_rhs(1,m,1,m),_uplo) );
635  assign( &nonconst_tri_ele(ex_gm_lhs(),_uplo), tri_ele(gms_rhs(1,m,1,m),_uplo) );
636  // Compute expected
637  for(int i = 1; i<=m; ++i) {
638  for(int j = i; j<=m; ++j) { // upper triangular
639  const int
640  _i = (_uplo == upper ? i : j), // adjust for lower triangular
641  _j = (_uplo == upper ? j : i);
642  ex_gm_lhs(_i,_j) += alpha * vs_rhs(i) * vs_rhs(j);
643  }
644  }
645  // Compute update
646  syr( alpha, vs_rhs(1,m), &nonconst_sym(gm_lhs(),_uplo) );
647  // Compare
648  update_success( result = comp( tri_ele(gm_lhs(),_uplo), tri_ele(ex_gm_lhs(),_uplo) )
649  , &success );
650  if( out && !result )
651  *out << "gm_lhs =\n" << gm_lhs << "ex_gm_lhs =\n" << ex_gm_lhs();
652  if(out) *out << " gm_lhs == expected_gm_lhs : " << result << std::endl;
653  }}
654 
655  // sym_lhs = alpha * vs_rhs1 * vs_rhs2' + alpha * vs_rhs2 * vs_rhs1' + sym_lhs (BLAS xSYR2).
656  if(out) *out << "\n****** sym_lhs = alpha * vs_rhs1 * vs_rhs2' + alpha * vs_rhs2 * vs_rhs1' + sym_lhs (BLAS xSYR2)\n";
657 
658  gm_lhs = 0.0;
659  ex_gm_lhs = 0.0;
660  alpha = 2.0;
661  DVector vs_rhs2 = vs_rhs.rev();
662 
663  {for( int i_uplo = 0; i_uplo < 2; ++i_uplo) {
664  const BLAS_Cpp::Uplo _uplo = a_uplo[i_uplo];
665  if(out)
666  *out << "\nFor sym_lhs = sym(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
667  << ")\n";
668  assign( &nonconst_tri_ele(gm_lhs(),_uplo), tri_ele(gms_rhs(1,m,1,m),_uplo) );
669  assign( &nonconst_tri_ele(ex_gm_lhs(),_uplo), tri_ele(gms_rhs(1,m,1,m),_uplo) );
670  // Compute expected
671  for(int i = 1; i<=m; ++i) {
672  for(int j = i; j<=m; ++j) { // upper triangular
673  const int
674  _i = (_uplo == upper ? i : j), // adjust for lower triangular
675  _j = (_uplo == upper ? j : i);
676  ex_gm_lhs(_i,_j) += alpha * (vs_rhs(i) * vs_rhs2(j) + vs_rhs2(i) * vs_rhs(j));
677  }
678  }
679  // Compute update
680  syr2( alpha, vs_rhs(1,m), vs_rhs2(1,m), &nonconst_sym(gm_lhs(),_uplo) );
681  // Compare
682  update_success( result = comp( tri_ele(gm_lhs(),_uplo), tri_ele(ex_gm_lhs(),_uplo) )
683  , &success );
684  if( out && !result )
685  *out << "gm_lhs =\n" << gm_lhs() << "ex_gm_lhs =\n" << ex_gm_lhs();
686  if(out) *out << " gm_lhs == expected_gm_lhs : " << result << std::endl;
687  }}
688 
689  // /////////////////////////////////////////
690  // Level-3 BLAS
691 
692  if(out) *out << "\n***\n*** Level-3 BLAS\n***\n";
693 
694  // Triangular Matrices
695  if(out) *out << "\n*** BLAS-2 Equivalent Triangular Matrix operations\n";
696 
697  DMatrix Tmp1, Tmp2;
698 
699  {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
700  for(int i_diag = 0; i_diag < 2; ++i_diag) {
701  for(int i_trans1 = 0; i_trans1 < 2; ++i_trans1) {
702  for(int i_trans2 = 0; i_trans2 < 2; ++i_trans2) {
703  const BLAS_Cpp::Transp
704  _trans1 = a_trans[i_trans1],
705  _trans2 = a_trans[i_trans2];
706  if(out)
707  *out << "\nLet M = tri(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
708  << "," << str_diag[i_diag] << ")"
709  << trans_char(_trans1) << std::endl;
710  DMatrixSliceTri M = tri(gms_rhs(1,m,1,m),a_uplo[i_uplo],a_diag[i_diag]);
711  // gm_lhs (left)
712  // gm_lhs = alpha * op(tri_rhs1) * op(gms_rhs2) (left) (BLAS xTRMM).
713  // gm_lhs = alpha * inv(op(tri_rhs1)) * op(gms_rhs2) (left) (BLAS xTRSM).
714  Tmp1.resize(1,1);
715  Tmp2.resize(1,1);
716  M_StInvMtM( &Tmp1, 2.0, M, _trans1, gms_rhs(1,m,1,m), _trans2 );
717  M_StMtM( &Tmp2, 0.5, M, _trans1, Tmp1(), no_trans );
718  update_success( result = comp( Tmp2(), no_trans, gms_rhs(1,m,1,m), _trans2 )
719  , &success );
720  if(out) *out << "(gm_lhs,left)...0.5*M*(2*inv(M)*gms_rhs(1,m,1,m)"
721  << trans_char(_trans2) << ")"
722  " == gms_rhs(1,m,1,m)" << trans_char(_trans2) << " : "
723  << result << std::endl;
724  // gms_lhs (left)
725  // gms_lhs = alpha * op(tri_rhs1) * op(gms_rhs2) (left) (BLAS xTRMM).
726  // gms_lhs = alpha * inv(op(tri_rhs1)) * op(gms_rhs2) (left) (BLAS xTRSM).
727  assign( &Tmp2, gms_rhs(1,m,1,m), _trans2 );
728  M_StInvMtM( &Tmp1(), 2.0, M, _trans1, Tmp2(), no_trans );
729  M_StMtM( &Tmp2, 0.5, M, _trans1, Tmp1(), no_trans );
730  update_success( result = comp( Tmp2(), no_trans, gms_rhs(1,m,1,m), _trans2 )
731  , &success );
732  if(out) {
733  *out
734  << "(gms_lhs,left)...0.5*M*(2*inv(M)*gms_rhs(1,m,1,m)"
735  << trans_char(_trans2) << ")"
736  " == gms_rhs(1,m,1,m)" << trans_char(_trans2) << " : "
737  << result << std::endl;
738  if(!result) {
739  *out
740  << "\ngms_lhs =\n" << Tmp2
741  << "\ngms_rhs(1,m,1,m) =\n" << gms_rhs(1,m,1,m);
742  }
743  }
744  // gm_lhs (right)
745  // gm_lhs = alpha * op(gms_rhs1) * op(tri_rhs2) (right) (BLAS xTRMM).
746  // gm_lhs = alpha * op(gms_rhs1) * inv(op(tri_rhs2)) (right) (BLAS xTRSM).
747  Tmp1.resize(1,1);
748  Tmp2.resize(1,1);
749  M_StMtInvM( &Tmp1, 2.0, gms_rhs(1,m,1,m), _trans1, M, _trans2 );
750  M_StMtM( &Tmp2, 0.5, Tmp1(), no_trans, M, _trans2 );
751  update_success( result = comp( Tmp2(), no_trans, gms_rhs(1,m,1,m), _trans1 )
752  , &success );
753  if(out) *out << "(gm_lhs,right)...5.0*(2.0*gms_rhs(1,m,1,m)"
754  << trans_char(_trans1) << "*inv(M))*M"
755  " == gms_rhs(1,m,1,m)" << trans_char(_trans1) << " : "
756  << result << std::endl;
757  // gms_lhs (right)
758  // gms_lhs = alpha * op(gms_rhs1) * op(tri_rhs2) (right) (BLAS xTRMM).
759  // gms_lhs = alpha * op(gms_rhs1) * inv(op(tri_rhs2)) (right) (BLAS xTRSM).
760  assign( &Tmp2(), gms_rhs(1,m,1,m), _trans1 );
761  M_StMtInvM( &Tmp2(), 2.0, Tmp2(), no_trans, M, _trans2 );
762  M_StMtM( &Tmp2(), 0.5, Tmp2(), no_trans, M, _trans2 );
763  update_success( result = comp( Tmp2(), no_trans, gms_rhs(1,m,1,m), _trans1 )
764  , &success );
765  if(out) *out << "(gms_lhs,right)...5.0*(2.0*gms_rhs(1,m,1,m)"
766  << trans_char(_trans1) << "*inv(M))*M"
767  " == gms_rhs(1,m,1,m)" << trans_char(_trans1) << " : "
768  << result << std::endl;
769  }
770  }
771  }
772  }}
773 
774  // *** LinAlgOpPack complient
775  if(out) *out << "\n*** Test DenseLinAlgPack complient (Mp_StMtM(...))\n";
776 
777  DMatrix I(0.0,n,n);
778  I.diag() = 1.0;
779 
780  // ****** Rectangular Matrices
781  if(out) *out << "\n****** Rectangular MtM (BLAS xGEMV)\n";
782 
783  // gms_lhs = alpha * op(gms_rhs1) * op(gms_rhs2) + beta * gms_lhs (BLAS xGEMV).
784  {for(int i_trans_B = 0; i_trans_B < 2; ++i_trans_B) {
785  for(int i_trans_I = 0; i_trans_I < 2; ++i_trans_I) {
786  const BLAS_Cpp::Transp
787  _trans_B = a_trans[i_trans_B],
788  _trans_I = a_trans[i_trans_I];
789  if(out) *out << "\nLet B = gms_rhs\n";
790  test_MtM( gms_rhs, _trans_B, I, _trans_I, out, &Tmp1, &Tmp2, &success );
791  }
792  }}
793 
794  // ****** Symmetric Matrices
795  if(out) *out << "\n****** Symmetric MtM\n";
796 
797  // gms_lhs = alpha * op(sym_rhs1) * op(gms_rhs2) + beta * gms_lhs (left) (BLAS xSYMM).
798  // gms_lhs = alpha * op(gms_rhs1) * op(sym_rhs2) + beta * gms_lhs (right) (BLAS xSYMM).
799  {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
800  for(int i_trans_B = 0; i_trans_B < 2; ++i_trans_B) {
801  for(int i_trans_I = 0; i_trans_I < 2; ++i_trans_I) {
802  const BLAS_Cpp::Uplo
803  _uplo = a_uplo[i_uplo];
804  const BLAS_Cpp::Transp
805  _trans_B = a_trans[i_trans_B],
806  _trans_I = a_trans[i_trans_I];
807  if(out) *out << "\nLet B = sym(gms_rhs(1,m,1,m)," << str_uplo[i_uplo] << ")\n";
808  test_MtM( sym(gms_rhs(1,m,1,m),_uplo), _trans_B, I, _trans_I, out
809  , &Tmp1, &Tmp2, &success );
810  }
811  }
812  }}
813 
814  // ****** Triangular Matrices
815  if(out) *out << "\n****** Triangular MtM\n";
816 
817  // gms_lhs = alpha * op(tri_rhs1) * op(gms_rhs2) + beta * gms_lhs (left).
818  // gms_lhs = alpha * op(gms_rhs1) * op(tri_rhs2) + beta * gms_lhs (right).
819 
820  {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
821  for(int i_diag = 0; i_diag < 2; ++i_diag) {
822  for(int i_trans_B = 0; i_trans_B < 2; ++i_trans_B) {
823  for(int i_trans_I = 0; i_trans_I < 2; ++i_trans_I) {
824  const BLAS_Cpp::Uplo
825  _uplo = a_uplo[i_uplo];
826  const BLAS_Cpp::Diag
827  _diag = a_diag[i_diag];
828  const BLAS_Cpp::Transp
829  _trans_B = a_trans[i_trans_B],
830  _trans_I = a_trans[i_trans_I];
831  if(out) *out << "\nLet B = tri(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
832  << "," << str_diag[i_diag] << ")\n";
833  test_MtM( tri(gms_rhs(1,m,1,m),_uplo,_diag), _trans_B, I, _trans_I, out
834  , &Tmp1, &Tmp2, &success );
835  }
836  }
837  }
838  }}
839 
840  // *** Symmetric Matrix updating
841  if(out) *out << "\n*** Symmetric Matrix updating\n";
842 
843  if(out) *out << "\nWarning! Not Tested!\n";
844 
845  // sym_lhs = alpha * op(gms_rhs) * op(gms_rhs') + beta * sym_lhs (BLAS xSYRK).
846  // sym_lhs = alpha * op(gms_rhs1) * op(gms_rhs2') + alpha * op(gms_rhs2) * op(gms_rhs1') + beta * sym_lhs (BLAS xSYR2K)
847 
848  } // end try
849  catch( const std::exception& excpt ) {
850  success = false;
851  if(out)
852  (*out) << "\nError, a standard exception was thrown: "
853  << typeName(excpt) << ": "
854  << excpt.what() << std::endl;
855  }
856  catch(...) {
857  success = false;
858  if(out)
859  (*out) << "\nError, an unknown exception was thrown\n";
860  }
861 
862  if(out) {
863  if(success)
864  (*out)
865  << "\n*** Congradulations, GenMatrixOp operations seem to check out. ***\n";
866  else
867  (*out)
868  << "\n*** Oops, all of the tests for GenMatrixOp operations "
869  "where not successful. ***\n";
870  }
871 
872 
873  return success;
874 }
LinAlgOpPack::M_StMtM
void M_StMtM(MatrixOp *M_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2)
M_lhs = alpha * op(M_rhs1) * op(M_rhs2).
Definition: AbstractLinAlgPack_AbstractLinAlgOpPack.cpp:103
DenseLinAlgPack::size_type
Teuchos::Ordinal size_type
Typedef for the size type of elements that are used by the library.
Definition: DenseLinAlgPack_Options.hpp:104
LinAlgOpPack::size_type
AbstractLinAlgPack::size_type size_type
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:54
typeName
std::string typeName(const T &t)
BLAS_Cpp::Uplo
Uplo
UPLO.
Definition: BLAS_Cpp_Types.hpp:79
DenseLinAlgPack::Vp_StV
void Vp_StV(DVectorSlice *vs_lhs, value_type alpha, const DVectorSlice &vs_rhs)
vs_lhs += alpha * vs_rhs (BLAS xAXPY)
Definition: DenseLinAlgPack_DVectorOpBLAS.cpp:120
DenseLinAlgPack::assign
void assign(DMatrix *gm_lhs, value_type alpha)
gm_lhs = alpha (elementwise)
Definition: DenseLinAlgPack_DMatrixOpBLAS.cpp:121
AbstractLinAlgPack::Mp_StMtM
void Mp_StMtM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2, value_type beta=1.0)
mwo_lhs = alpha * op(mwo_rhs1) * op(mwo_rhs2) + beta * mwo_lhs (right) (xGEMM).
Definition: AbstractLinAlgPack_MatrixOp.cpp:579
AbstractLinAlgPack::Vt_S
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
Definition: AbstractLinAlgPack_VectorStdOps.cpp:176
DenseLinAlgPack::DMatrix::diag
DVectorSlice diag(difference_type k=0)
Definition: DenseLinAlgPack_DMatrixClass.hpp:1082
BLAS_Cpp::trans_to_bool
bool trans_to_bool(Transp _trans)
Returns true if _trans == trans.
Definition: BLAS_Cpp_Types.hpp:97
DenseLinAlgPack::tri_ele
const DMatrixSliceTriEle tri_ele(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo)
Definition: DenseLinAlgPack_DMatrixAsTriSym.hpp:119
AbstractLinAlgPack::Vp_StV
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
Definition: AbstractLinAlgPack_VectorStdOps.cpp:191
BLAS_Cpp::rows
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return rows of a possible transposed matrix.
Definition: BLAS_Cpp_Types.hpp:120
BLAS_Cpp::syr2
void syr2(Uplo uplo, f_int n, f_dbl_prec alpha, const f_dbl_prec *px, f_int incx, const f_dbl_prec *py, f_int incy, f_dbl_prec *pa, f_int lda)
Definition: DenseLinAlgPack_BLAS_Cpp.cpp:471
AbstractLinAlgPack::M_StMtInvM
void M_StMtInvM(MatrixOp *m_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixNonsing &M_rhs2, BLAS_Cpp::Transp trans_rhs2)
m_lhs = alpha * op(mwo_rhs1) * inv(op(M_rhs2)) (left)
Definition: AbstractLinAlgPack_MatrixNonsing.hpp:296
LinAlgOpPack::assign
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:442
DenseLinAlgPack::DVectorSlice
VectorSliceTmpl< value_type > DVectorSlice
Definition: DenseLinAlgPack_Types.hpp:94
Teuchos::SerialDenseVector::resize
int resize(OrdinalType length_in)
DenseLinAlgPack::VectorTmpl::resize
void resize(size_type n, value_type val=value_type())
Definition: DenseLinAlgPack_DVectorClassTmpl.hpp:1095
DenseLinAlgPack_MatVecCompare.hpp
DenseLinAlgPack_DMatrixOut.hpp
DenseLinAlgPack::sym
const DMatrixSliceSym sym(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo)
Definition: DenseLinAlgPack_DMatrixAsTriSym.hpp:258
BLAS_Cpp::trans
Transposed.
Definition: BLAS_Cpp_Types.hpp:75
DenseLinAlgPack::tri
const DMatrixSliceTri tri(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo, BLAS_Cpp::Diag diag)
Definition: DenseLinAlgPack_DMatrixAsTriSym.hpp:192
DenseLinAlgPack_DMatrixOp.hpp
BLAS_Cpp::lower
Definition: BLAS_Cpp_Types.hpp:79
BLAS_Cpp::no_trans
Not transposed.
Definition: BLAS_Cpp_Types.hpp:74
DenseLinAlgPack_TestDenseLinAlgPack.hpp
BLAS_Cpp::unit
Definition: BLAS_Cpp_Types.hpp:84
DenseLinAlgPack::nonconst_sym
DMatrixSliceSym nonconst_sym(DMatrixSlice gms, BLAS_Cpp::Uplo uplo)
Return a symmetric matrix.
Definition: DenseLinAlgPack_DMatrixAsTriSym.hpp:251
QPKWIK_Output::out
std::ostream * out
BLAS_Cpp::upper
Definition: BLAS_Cpp_Types.hpp:79
M
const LAPACK_C_Decl::f_int & M
Definition: DenseLinAlgPack_LAPACK_C_Decl.cpp:82
DenseLinAlgPack::Mp_StMtM
void Mp_StMtM(DMatrixSlice *gms_lhs, value_type alpha, const DMatrixSlice &gms_rhs1, BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice &gms_rhs2, BLAS_Cpp::Transp trans_rhs2, value_type beta=1.0)
Definition: DenseLinAlgPack_DMatrixOpBLAS.cpp:476
AbstractLinAlgPack::V_InvMtV
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
v_lhs = inv(op(M_rhs1)) * v_rhs2
Definition: AbstractLinAlgPack_MatrixNonsing.hpp:249
AbstractLinAlgPack::Mp_StM
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
Definition: AbstractLinAlgPack_MatrixOp.cpp:446
BLAS_Cpp::Diag
Diag
DIAG.
Definition: BLAS_Cpp_Types.hpp:84
AbstractLinAlgPack::Vp_StMtV
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
Definition: AbstractLinAlgPack_MatrixOp.hpp:857
DenseLinAlgPack_DMatrixClass.hpp
DenseLinAlgPack::DMatrix::resize
void resize(size_type rows, size_type cols, value_type val=value_type())
Resize matrix to a (rows x cols) matrix and initializes any added elements by val.
Definition: DenseLinAlgPack_DMatrixClass.hpp:1010
DenseLinAlgPack::nonconst_tri_ele
DMatrixSliceTriEle nonconst_tri_ele(DMatrixSlice gms, BLAS_Cpp::Uplo uplo)
Return a triangular element-wise matrix.
Definition: DenseLinAlgPack_DMatrixAsTriSym.hpp:112
AbstractLinAlgPack::Mt_S
void Mt_S(MatrixOp *mwo_lhs, value_type alpha)
Definition: AbstractLinAlgPack_MatrixOp.cpp:438
BLAS_Cpp::syr
void syr(Uplo uplo, f_int n, f_dbl_prec alpha, const f_dbl_prec *px, f_int incx, f_dbl_prec *pa, f_int lda)
Definition: DenseLinAlgPack_BLAS_Cpp.cpp:455
AbstractLinAlgPack::DVectorSlice
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
Definition: AbstractLinAlgPack_Types.hpp:67
DenseLinAlgPack::Vt_S
void Vt_S(DVectorSlice *vs_lhs, value_type alpha)
vs_lhs *= alpha (BLAS xSCAL) (*** Note that alpha == 0.0 is handeled as vs_lhs = 0.0)
Definition: DenseLinAlgPack_DVectorOpBLAS.cpp:111
LinAlgOpPack::V_MtV
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:500
LinAlgOpPack::value_type
AbstractLinAlgPack::value_type value_type
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:56
DenseLinAlgPack::DMatrixSlice::diag
DVectorSlice diag(difference_type k=0)
Definition: DenseLinAlgPack_DMatrixClass.hpp:878
DenseLinAlgPack::Mp_StM
void Mp_StM(DMatrixSliceTriEle *tri_lhs, value_type alpha, const DMatrixSliceTriEle &tri_rhs)
tri_lhs += alpha * tri_rhs (BLAS xAXPY)
Definition: DenseLinAlgPack_DMatrixOpBLAS.cpp:245
DenseLinAlgPack::value_type
FortranTypes::f_dbl_prec value_type
Typedef for the value type of elements that is used for the library.
Definition: DenseLinAlgPack_Options.hpp:100
DenseLinAlgPack::DVector
VectorTmpl< value_type > DVector
Definition: DenseLinAlgPack_Types.hpp:90
BLAS_Cpp::nonunit
Definition: BLAS_Cpp_Types.hpp:84
DenseLinAlgPack::Vp_StMtV
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, const DMatrixSlice &gms_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta=1.0)
vs_lhs = alpha * op(gms_rhs1) * vs_rhs2 + beta * vs_lhs (BLAS xGEMV)
Definition: DenseLinAlgPack_DMatrixOpBLAS.cpp:349
DenseLinAlgPack_DVectorOut.hpp
TestingHelperPack::update_success
bool update_success(bool result_check, bool *success)
Helper function for updating a flag for if an operation returned false.
Definition: TestingHelperPack_update_success.cpp:48
AbstractLinAlgPack::M_StInvMtM
void M_StInvMtM(MatrixOp *m_lhs, value_type alpha, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2)
m_lhs = alpha * inv(op(mwo_rhs1)) * op(mwo_rhs2) (right)
Definition: AbstractLinAlgPack_MatrixNonsing.hpp:287
DenseLinAlgPack::comp
bool comp(const DVectorSlice &vs1, const DVectorSlice &vs2)
Definition: DenseLinAlgPack_MatVecCompare.cpp:72
AbstractLinAlgPack::DMatrixSlice
DenseLinAlgPack::DMatrixSlice DMatrixSlice
Definition: AbstractLinAlgPack_Types.hpp:71
BLAS_Cpp::Transp
Transp
TRANS.
Definition: BLAS_Cpp_Types.hpp:73
AbstractLinAlgPack::DMatrix
DenseLinAlgPack::DMatrix DMatrix
Definition: AbstractLinAlgPack_Types.hpp:69
n
int n
BLAS_Cpp::cols
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return columns of a possible transposed matrix.
Definition: BLAS_Cpp_Types.hpp:127
Generated on Wed Dec 2 2015 08:50:25 for MOOCHO (Single Doxygen Collection) by   doxygen 1.8.6