Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_OskiMatrix.cpp
Go to the documentation of this file.
1
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42
43#include "Epetra_ConfigDefs.h"
44#include "Epetra_Map.h"
45
46#ifdef HAVE_OSKI
47#ifdef HAVE_EPETRA_TEUCHOS
48#include "Epetra_OskiMatrix.h"
49#include "Epetra_Import.h"
50
51//=============================================================================
52
54 : Epetra_CrsMatrix(Source),
55 Epetra_View_(Source.Epetra_View_),
56 Copy_Created_(true) {
57 A_tunable_ = oski_CopyMat(Source.A_tunable_);
58}
59
61 const Teuchos::ParameterList& List)
62 : Epetra_CrsMatrix(Source),
63 Epetra_View_(&Source) {
64 bool AutoTune = false;
65 bool DeepCopy = false;
66 string Matrix = "general\0";
67 bool IsDiagNotStored = false;
68 bool IsArrayZeroBased = false;
69 int MyIndexBase = 1;
70 bool AreIndicesSorted = false;
71 bool AreIndicesRepeated = false;
72 oski_inmatprop_t MatrixType = MAT_GENERAL;
73 oski_inmatprop_t Diagonal = MAT_DIAG_EXPLICIT;
74 oski_inmatprop_t ArrayBasis = INDEX_ONE_BASED;
75 oski_inmatprop_t SortedIndices = INDEX_UNSORTED;
76 oski_inmatprop_t RepeatedIndices = INDEX_REPEATED;
77 int* RowPtr = NULL;
78 int* IndPtr = NULL;
79 double* ValPtr = NULL;
80 AutoTune = const_cast <Teuchos::ParameterList &>(List).get("autotune", false);
81 if(List.isParameter("deepcopy"))
82 DeepCopy = Teuchos::getParameter<bool>(List, "deepcopy");
83 if(AutoTune){ //Use parameters from the Epetra matrix to set as many fields as possible
84 if(LowerTriangular())
85 MatrixType = MAT_TRI_LOWER;
86 if(UpperTriangular())
87 MatrixType = MAT_TRI_UPPER;
88 if(Sorted())
89 SortedIndices = INDEX_SORTED;
90 MyIndexBase = IndexBase();
91 if(MyIndexBase == 0)
92 ArrayBasis = INDEX_ZERO_BASED;
93 else if(MyIndexBase == 1)
94 ArrayBasis = INDEX_ONE_BASED;
95 else if(!List.isParameter("zerobased")) {
96 std::cerr << "An OskiMatrix must be either zero or one based.\n";
97 return;
98 }
99 if(NoRedundancies())
100 RepeatedIndices = INDEX_UNIQUE;
101 }
102 if(List.isParameter("matrixtype")) {
103 Matrix = Teuchos::getParameter<string>(List, "matrixtype");
104 if(Matrix == "general")
105 MatrixType = MAT_GENERAL;
106 else if(Matrix == "uppertri")
107 MatrixType = MAT_TRI_UPPER;
108 else if(Matrix == "lowertri")
109 MatrixType = MAT_TRI_LOWER;
110 else if(Matrix == "uppersymm")
111 MatrixType = MAT_SYMM_UPPER;
112 else if(Matrix == "lowersymm")
113 MatrixType = MAT_SYMM_LOWER;
114 else if(Matrix == "fullsymm")
115 MatrixType = MAT_SYMM_FULL;
116 else if(Matrix == "upperherm")
117 MatrixType = MAT_HERM_UPPER;
118 else if(Matrix == "lowerherm")
119 MatrixType = MAT_HERM_LOWER;
120 else if(Matrix == "fullherm")
121 MatrixType = MAT_HERM_FULL;
122 }
123 if(List.isParameter("diagstored"))
124 IsDiagNotStored = Teuchos::getParameter<bool>(List, "diagstored");
125 if(List.isParameter("zerobased"))
126 IsArrayZeroBased = Teuchos::getParameter<bool>(List, "zerobased");
127 if(List.isParameter("sorted"))
128 AreIndicesSorted = Teuchos::getParameter<bool>(List, "sorted");
129 if(List.isParameter("unique"))
130 DeepCopy = Teuchos::getParameter<bool>(List, "unique");
131 if(IsDiagNotStored)
132 Diagonal = MAT_UNIT_DIAG_IMPLICIT;
133 if(IsArrayZeroBased)
134 ArrayBasis = INDEX_ZERO_BASED;
135 if(AreIndicesSorted)
136 SortedIndices = INDEX_SORTED;
137 if(AreIndicesRepeated)
138 RepeatedIndices = INDEX_UNIQUE;
139 if(ExtractCrsDataPointers(RowPtr, IndPtr, ValPtr)) {
140 std::cerr << "Cannot wrap matrix as an Oski Matrix because at least one of FillComplete and Optimize Storage has not been called\n";
141 return;
142 }
143 if(DeepCopy) {
144 Copy_Created_ = true;
145 A_tunable_ = oski_CreateMatCSR(RowPtr, IndPtr, ValPtr, Source.NumMyRows(), Source.NumMyCols(), COPY_INPUTMAT, 5, MatrixType, Diagonal, ArrayBasis, SortedIndices, RepeatedIndices);
146 }
147 else {
148 Copy_Created_ = false;
149 A_tunable_ = oski_CreateMatCSR(RowPtr, IndPtr, ValPtr, Source.NumMyRows(), Source.NumMyCols(), SHARE_INPUTMAT, 5, MatrixType, Diagonal, ArrayBasis, SortedIndices, RepeatedIndices);
150 }
151}
152
154 if(oski_DestroyMat(A_tunable_))
155 std::cerr << "Destroy Matrix failed.\n";
156}
157
159 int NumEntries,
160 double* Values,
161 int* Indices) {
162 int ReturnVal = 0;
163 if (Copy_Created_) {
164 for(int i = 0; i < NumEntries; i++) {
165 ReturnVal = oski_SetMatEntry(A_tunable_, MyRow, Indices[i], Values[i]);
166 if(ReturnVal)
167 break;
168 }
169 if(!ReturnVal)
170 ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->ReplaceMyValues(MyRow, NumEntries, Values, Indices);
171 }
172 else
173 ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->ReplaceMyValues(MyRow, NumEntries, Values, Indices);
174 if(ReturnVal)
175 std::cerr << "Error in ReplaceMyValues\n";
176 return ReturnVal;
177}
178
180 int NumEntries,
181 double* Values,
182 int* Indices) {
183 int ReturnVal = 0;
184 if (Copy_Created_) {
185 for(int i = 0; i < NumEntries; i++) {
186 ReturnVal = oski_SetMatEntry(A_tunable_, MyRow, Indices[i], oski_GetMatEntry(A_tunable_, MyRow, Indices[i]));
187 if(ReturnVal)
188 break;
189 }
190 if(!ReturnVal)
191 ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->SumIntoMyValues(MyRow, NumEntries, Values, Indices);
192 }
193 else
194 ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->SumIntoMyValues(MyRow, NumEntries, Values, Indices);
195 if(ReturnVal)
196 std::cerr << "Error in SumIntoMyValues\n";
197 return ReturnVal;
198}
199
201 int ReturnValue = 0;
202 ReturnValue = Epetra_View_->ExtractDiagonalCopy(Diagonal);
203 if (ReturnValue)
204 std::cerr << "Error in ExtractDiagonalCopy\n";
205 return ReturnValue;
206}
207
209 int ReturnVal = 0;
210 if (Copy_Created_) {
211 ReturnVal = oski_SetMatDiagValues(A_tunable_, 0, Diagonal.Oski_View());
212 if(!ReturnVal)
213 ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->ReplaceDiagonalValues(*Diagonal.Epetra_View());
214 }
215 else
216 ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->ReplaceDiagonalValues(*Diagonal.Epetra_View());
217 if(ReturnVal)
218 std::cerr << "Error in ReplaceDiagonalValues\n";
219 return ReturnVal;
220}
221
222int Epetra_OskiMatrix::Multiply(bool TransA,
223 const Epetra_Vector& x,
224 Epetra_Vector& y) const {
225 int ReturnVal;
226 ReturnVal = this->Multiply(TransA, x, y, 1.0, 0.0);
227 return ReturnVal;
228}
229
230int Epetra_OskiMatrix::Multiply(bool TransA,
231 const Epetra_Vector& x,
232 Epetra_Vector& y,
233 double Alpha,
234 double Beta) const {
235 int ReturnVal;
236
237 if(!Filled())
238 EPETRA_CHK_ERR(-1); // Matrix must be filled.
239
240 double* xp = (double*) x.Values();
241 double* yp = (double*) y.Values();
242
243 Epetra_Vector * xcopy = 0;
244 if (&x==&y && Importer()==0 && Exporter()==0) {
245 xcopy = new Epetra_Vector(x);
246 xp = (double *) xcopy->Values();
247 }
248 UpdateImportVector(1); // Refresh import and output vectors if needed
250
251 if(!TransA) {
252
253 // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
254 if(Importer() != 0) {
256 xp = (double*) ImportVector_->Values();
257 }
258
259 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
260 if(Exporter() != 0)
261 yp = (double*) ExportVector_->Values();
262
263
264 oski_vecview_t oskiX;
265 oski_vecview_t oskiY;
266 if(Importer() != 0)
267 oskiX = oski_CreateVecView(xp,ImportVector_->MyLength(),1);
268 else
269 oskiX = oski_CreateVecView(xp,x.MyLength(),1);
270 if(Exporter() != 0)
271 oskiY = oski_CreateVecView(yp,ExportVector_->MyLength(),1);
272 else
273 oskiY = oski_CreateVecView(yp,y.MyLength(),1);
274
275 //Do actual computation
276 ReturnVal = oski_MatMult(A_tunable_, OP_NORMAL, Alpha, oskiX, Beta, oskiY);
277
278 if(Exporter() != 0) {
279 y.PutScalar(0.0); // Make sure target is zero
280 EPETRA_CHK_ERR(y.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
281 }
282 // Handle case of rangemap being a local replicated map
283 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
284 } //if(!TransA)
285 else { // Transpose operation
286
287 // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
288 if(Exporter() != 0) {
290 xp = (double*) ExportVector_->Values();
291 }
292
293 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
294 if(Importer() != 0)
295 yp = (double*) ImportVector_->Values();
296
297 oski_vecview_t oskiX;
298 oski_vecview_t oskiY;
299 if(Importer() != 0)
300 oskiY = oski_CreateVecView(yp,ImportVector_->MyLength(),1);
301 else
302 oskiY = oski_CreateVecView(yp,y.MyLength(),1);
303 if(Exporter() != 0)
304 oskiX = oski_CreateVecView(xp,ExportVector_->MyLength(),1);
305 else
306 oskiX = oski_CreateVecView(xp,x.MyLength(),1);
307
308 // Do actual computation
309 ReturnVal = oski_MatMult(A_tunable_, OP_TRANS, Alpha, oskiX, Beta, oskiY);
310
311 if(Importer() != 0) {
312 y.PutScalar(0.0); // Make sure target is zero
313 EPETRA_CHK_ERR(y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
314 }
315 // Handle case of rangemap being a local replicated map
316 if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
317 }
318 if(ReturnVal)
319 std::cerr << "OskiVector multiply error\n";
321 return ReturnVal;
322}
323
324int Epetra_OskiMatrix::Multiply(bool TransA,
325 const Epetra_MultiVector& X,
326 Epetra_MultiVector& Y) const {
327 int ReturnVal;
328 ReturnVal = this->Multiply(TransA, X, Y, 1.0, 0.0);
329 return ReturnVal;
330}
331
332int Epetra_OskiMatrix::Multiply(bool TransA,
333 const Epetra_MultiVector& X,
335 double Alpha,
336 double Beta) const {
337 int ReturnVal;
338
339 if(!Filled())
340 EPETRA_CHK_ERR(-1); // Matrix must be filled.
341
342 int NumVectors = X.NumVectors();
343 if (NumVectors!=Y.NumVectors()) {
344 EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
345 }
346
347 double** Xp = (double**) X.Pointers();
348 double** Yp = (double**) Y.Pointers();
349
350 int LDX = X.ConstantStride() ? X.Stride() : 0;
351 int LDY = Y.ConstantStride() ? Y.Stride() : 0;
352
353 Epetra_MultiVector* Xcopy = 0;
354 if (&X==&Y && Importer()==0 && Exporter()==0) {
355 Xcopy = new Epetra_MultiVector(X);
356 Xp = (double **) Xcopy->Pointers();
357 LDX = Xcopy->ConstantStride() ? Xcopy->Stride() : 0;
358 }
359 UpdateImportVector(NumVectors); // Refresh import and output vectors if needed
360 UpdateExportVector(NumVectors);
361
362 if (!TransA) {
363
364 // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
365 if (Importer()!=0) {
367 Xp = (double**)ImportVector_->Pointers();
369 }
370
371 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
372 if (Exporter()!=0) {
373 Yp = (double**)ExportVector_->Pointers();
375 }
376
377 oski_vecview_t oskiX;
378 oski_vecview_t oskiY;
379 if (Importer()!=0)
380 oskiX = oski_CreateMultiVecView(*Xp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
381 else
382 oskiX = oski_CreateMultiVecView(*Xp,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
383 if (Exporter()!=0)
384 oskiY = oski_CreateMultiVecView(*Yp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
385 else
386 oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
387 // Do actual computation
388 ReturnVal = oski_MatMult(A_tunable_, OP_NORMAL, Alpha, oskiX, Beta, oskiY);
389 if(ReturnVal)
390 std::cerr << "OskiMultiVector multiply error\n";
391 if (Exporter()!=0) {
392 Y.PutScalar(0.0); // Make sure target is zero
393 Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
394 }
395 // Handle case of rangemap being a local replicated map
396 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
397 }
398 else { // Transpose operation
399
400 // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
401
402 if (Exporter()!=0) {
404 Xp = (double**)ExportVector_->Pointers();
406 }
407
408 // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
409 if (Importer()!=0) {
410 Yp = (double**)ImportVector_->Pointers();
412 }
413
414 oski_vecview_t oskiX;
415 oski_vecview_t oskiY;
416 if (Importer()!=0)
417 oskiY = oski_CreateMultiVecView(*Yp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
418 else
419 oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
420 if (Exporter()!=0)
421 oskiX = oski_CreateMultiVecView(*Xp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
422 else
423 oskiX = oski_CreateMultiVecView(*Xp,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
424
425 // Do actual computation
426 ReturnVal = oski_MatMult(A_tunable_, OP_TRANS, Alpha, oskiX, Beta, oskiY);
427 if(ReturnVal)
428 std::cerr << "OskiMultiVector multiply error\n";
429 if (Importer()!=0) {
430 Y.PutScalar(0.0); // Make sure target is zero
431 EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
432 }
433 // Handle case of rangemap being a local replicated map
434 if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
435 }
437 //Y.ResetView(Yp);
438 return ReturnVal;
439}
440
441int Epetra_OskiMatrix::Solve(bool Upper, bool TransA, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const {
442 int ReturnVal;
443 ReturnVal = this->Solve(TransA, x, y, 1.0);
444 return ReturnVal;
445}
446
447int Epetra_OskiMatrix::Solve(bool TransA, const Epetra_Vector& x, Epetra_Vector& y, double Alpha) const {
448 std::cout << "This function Epetra_OskiMatrix::Solve probably works in serial but has not been tested.\n It will not work in parallel.\n If you wish to use it feel free to comment out this line and the next return statement.\n However, correctness and performance are not guaranteed.\n";
449 return(-1);
450 Epetra_OskiVector* xCast = NULL;
451 Epetra_OskiVector* yCast = NULL;
452 Epetra_OskiVector* tCast = NULL;
453 bool xCreate = false;
454 bool yCreate = false;
455 int ReturnVal;
456 xCast = dynamic_cast<Epetra_OskiVector*>(const_cast<Epetra_Vector*>(&x));
457 yCast = dynamic_cast<Epetra_OskiVector*>(&y);
458 if (xCast == NULL) {
459 xCast = new Epetra_OskiVector(x);
460 xCreate = true;
461 }
462 if (yCast == NULL) {
463 yCast = new Epetra_OskiVector(y);
464 yCreate = true;
465 }
466 tCast = new Epetra_OskiVector(x);
467 if(TransA)
468 ReturnVal = oski_MatTrisolve(A_tunable_, OP_TRANS, Alpha, (*tCast).Oski_View());
469 else
470 ReturnVal = oski_MatTrisolve(A_tunable_, OP_NORMAL, Alpha, (*tCast).Oski_View());
471 if(ReturnVal)
472 std::cerr << "OskiVector Solve error\n";
473 if(xCreate)
474 delete xCast;
475 yCast = tCast;
476 if(yCreate)
477 delete yCast;
478 delete tCast;
479 return ReturnVal;
480}
481
482int Epetra_OskiMatrix::Solve(bool Upper, bool TransA, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
483 int ReturnVal;
484 ReturnVal = this->Solve(TransA, X, Y, 1.0);
485 return ReturnVal;
486}
487
488int Epetra_OskiMatrix::Solve(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y, double Alpha) const {
489 std::cout << "This function Epetra_OskiMatrix::Solve probably works in serial but has not been tested.\n It will not work in parallel.\n If you wish to use it feel free to comment out this line and the next return statement.\n However, correctness and performance are not guaranteed.\n";
490 return(-1);
491 Epetra_OskiMultiVector* XCast = NULL;
492 Epetra_OskiMultiVector* YCast = NULL;
493 Epetra_OskiMultiVector* TCast = NULL;
494 bool XCreate = false;
495 bool YCreate = false;
496 int ReturnVal;
497 XCast = dynamic_cast<Epetra_OskiMultiVector*>(const_cast<Epetra_MultiVector*>(&X));
498 YCast = dynamic_cast<Epetra_OskiMultiVector*>(&Y);
499 if (XCast == NULL) {
500 XCast = new Epetra_OskiMultiVector(X);
501 XCreate = true;
502 }
503 if (YCast == NULL) {
504 YCast = new Epetra_OskiMultiVector(Y);
505 YCreate = true;
506 }
507 TCast = new Epetra_OskiMultiVector(X);
508 if(TransA)
509 ReturnVal = oski_MatTrisolve(A_tunable_, OP_TRANS, Alpha, (*TCast).Oski_View());
510 else
511 ReturnVal = oski_MatTrisolve(A_tunable_, OP_NORMAL, Alpha, (*TCast).Oski_View());
512 if(ReturnVal)
513 std::cerr << "OskiMultiVector Solve error\n";
514 if(XCreate)
515 delete XCast;
516 YCast = TCast;
517 if(YCreate)
518 delete YCast;
519 delete TCast;
520 return ReturnVal;
521}
522
524 const Epetra_Vector& x,
525 Epetra_Vector& y,
526 Epetra_Vector* t,
527 double Alpha,
528 double Beta) const {
529 int ReturnVal;
530
531 if(!Filled())
532 EPETRA_CHK_ERR(-1); // Matrix must be filled.
533
534 double* xp = (double*) x.Values();
535 double* xp2 = (double*) x.Values();
536 double* yp = (double*) y.Values();
537 double* tp = 0;
538 if(t != NULL)
539 tp = (double*) t->Values();
540
541 Epetra_Vector* xcopy = 0;
542 if (&x==&y && Importer()==0 && Exporter()==0) {
543 xcopy = new Epetra_Vector(x);
544 xp = (double *) xcopy->Values();
545 }
546 UpdateImportVector(1); // Refresh import and output vectors if needed
548
549
550 if(ATA) {
551 if(Importer() != 0) {
553 xp = (double*) ImportVector_->Values();
554 xp2 = new double[ImportVector_->MyLength()];
555 for(int i = 0; i < ImportVector_->MyLength(); i++)
556 xp2[i] = xp[i];
558 yp = (double*) ImportVector_->Values();
559 }
560
561 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
562 if(Exporter() != 0 && (t) != NULL) {
563 tp = (double*) ExportVector_->Values();
564 }
565
566 oski_vecview_t oskiX=0;
567 oski_vecview_t oskiY=0;
568 oski_vecview_t oskiT=0;
569 if(Importer() != 0)
570 oskiX = oski_CreateVecView(xp2,ImportVector_->MyLength(),1);
571 else
572 oskiX = oski_CreateVecView(xp2,x.MyLength(),1);
573 if(Importer() != 0)
574 oskiY = oski_CreateVecView(yp,ImportVector_->MyLength(),1);
575 else
576 oskiY = oski_CreateVecView(yp,y.MyLength(),1);
577
578 if (t != NULL) {
579 if(Exporter() != 0)
580 oskiT = oski_CreateVecView(tp,ExportVector_->MyLength(),1);
581 else
582 oskiT = oski_CreateVecView(tp,t->MyLength(),1);
583
584 }
585 else
586 oskiT = INVALID_VEC;
587 ReturnVal = oski_MatTransMatMult(A_tunable_, OP_AT_A, Alpha, oskiX, Beta, oskiY, oskiT);
588
589 if(Importer() != 0) {
590 y.PutScalar(0.0); // Make sure target is zero
591 EPETRA_CHK_ERR(y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
592 }
593 // Handle case of rangemap being a local replicated map
594 if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
595
596 if(Exporter() != 0 && (t != NULL)) {
597 t->PutScalar(0.0); // Make sure target is zero
598 EPETRA_CHK_ERR(t->Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
599 }
600 // Handle case of rangemap being a local replicated map
601 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(t->Reduce());
603
604 }
605 else {
606 if(this->Comm().NumProc() == 1) {
607 oski_vecview_t oskiX=0;
608 oski_vecview_t oskiY=0;
609 oski_vecview_t oskiT=0;
610 if (t != NULL)
611 oskiT = oski_CreateVecView(tp,t->MyLength(),1);
612 oskiX = oski_CreateVecView(xp,x.MyLength(),1);
613 oskiY = oski_CreateVecView(yp,y.MyLength(),1);
614 ReturnVal = oski_MatTransMatMult(A_tunable_, OP_A_AT, Alpha, oskiX, Beta, oskiY, oskiT);
616 }
617 else {
618
619 if(t == NULL) {
620 Epetra_Vector tempResult(this->DomainMap());
621 ReturnVal = this->Multiply(false, x, tempResult, 1.0, 0.0);
622 int ReturnVal2 = this->Multiply(true, tempResult, y, Alpha, Beta);
623 if(ReturnVal < ReturnVal2)
624 ReturnVal = ReturnVal2;
625 }
626 else {
627 ReturnVal = this->Multiply(false, x, *t, 1.0, 0.0);
628 int ReturnVal2 = this->Multiply(true,*t, y, Alpha, Beta);
629 if(ReturnVal < ReturnVal2)
630 ReturnVal = ReturnVal2;
631 }
632 }
633 }
634 if (xcopy!=0) {
635 delete xcopy;
636 EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of x
637 return(1);
638 }
639 return ReturnVal;
640}
641
643 const Epetra_MultiVector& X,
646 double Alpha,
647 double Beta) const {
648 int ReturnVal;
649
650 if(!Filled())
651 EPETRA_CHK_ERR(-1); // Matrix must be filled.
652
653 int NumVectors = X.NumVectors();
654 if (NumVectors!=Y.NumVectors()) {
655 EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
656 }
657
658
659
660 double** Xp = (double**) X.Pointers();
661 double** Xp2 = (double**) X.Pointers();
662 double** Yp = (double**) Y.Pointers();
663 double** Tp = 0;
664 if(T != NULL)
665 Tp = (double**) T->Pointers();
666
667 int LDX = X.ConstantStride() ? X.Stride() : 0;
668 int LDY = Y.ConstantStride() ? Y.Stride() : 0;
669 int LDT = 0;
670 if(T != NULL)
671 LDT = T->ConstantStride() ? T->Stride() : 0;
672
673 Epetra_MultiVector* Xcopy = 0;
674 Epetra_MultiVector* X2 = 0;
675 if (&X==&Y && Importer()==0 && Exporter()==0) {
676 Xcopy = new Epetra_MultiVector(X);
677 Xp = (double **) Xcopy->Pointers();
678 LDX = Xcopy->ConstantStride() ? Xcopy->Stride() : 0;
679 }
680 UpdateImportVector(NumVectors); // Refresh import and output vectors if needed
681 UpdateExportVector(NumVectors);
682
683
684 if(ATA) {
685 if(Importer() != 0) {
687 Xp = (double**) ImportVector_->Pointers();
689// need to think about this
690 X2 = new Epetra_MultiVector(X);
691 double** Xp2 = (double**) X2->Pointers();
692 Xp2 = (double **) X2->Pointers();
694 Yp = (double**) ImportVector_->Pointers();
696 }
697
698 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
699 if(Exporter() != 0 && T != NULL) {
700 Tp = (double**) ExportVector_->Pointers();
702 }
703 oski_vecview_t oskiX=0;
704 oski_vecview_t oskiY=0;
705 oski_vecview_t oskiT=0;
706 if(Importer() != 0)
707 oskiX = oski_CreateMultiVecView(*Xp2,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
708 else
709 oskiX = oski_CreateMultiVecView(*Xp2,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
710 if(Importer() != 0)
711 oskiY = oski_CreateMultiVecView(*Yp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
712 else
713 oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
714
715 if (T != NULL) {
716 if(Exporter() != 0)
717 oskiT = oski_CreateMultiVecView(*Tp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDT);
718 else
719 oskiT = oski_CreateMultiVecView(*Tp,T->MyLength(),NumVectors,LAYOUT_COLMAJ,LDT);
720
721 }
722 else
723 oskiT = INVALID_VEC;
724 ReturnVal = oski_MatTransMatMult(A_tunable_, OP_AT_A, Alpha, oskiX, Beta, oskiY, oskiT);
725
726 if(Importer() != 0) {
727 Y.PutScalar(0.0); // Make sure target is zero
728 EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
729 }
730 // Handle case of rangemap being a local replicated map
731 if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
732
733 if(Exporter() != 0 && (T != NULL)) {
734 T->PutScalar(0.0); // Make sure target is zero
735 EPETRA_CHK_ERR(T->Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
736 }
737 // Handle case of rangemap being a local replicated map
738 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(T->Reduce());
740
741 }
742 else {
743 if(this->Comm().NumProc() == 1) {
744 oski_vecview_t oskiX=0;
745 oski_vecview_t oskiY=0;
746 oski_vecview_t oskiT=0;
747 if (T != NULL)
748 oskiT = oski_CreateMultiVecView(*Tp,T->MyLength(),NumVectors, LAYOUT_COLMAJ,LDT);
749 oskiX = oski_CreateMultiVecView(*Xp,X.MyLength(),NumVectors, LAYOUT_COLMAJ,LDX);
750 oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors, LAYOUT_COLMAJ,LDY);
751 ReturnVal = oski_MatTransMatMult(A_tunable_, OP_A_AT, Alpha, oskiX, Beta, oskiY, oskiT);
752 UpdateFlops(4 * NumGlobalNonzeros64() *NumVectors);
753 }
754 else {
755 if(T == NULL) {
756 Epetra_MultiVector TempResult(this->DomainMap(), NumVectors);
757 ReturnVal = this->Multiply(false, X, TempResult, 1.0, 0.0);
758 int ReturnVal2 = this->Multiply(true, TempResult, Y, Alpha, Beta);
759 if(ReturnVal < ReturnVal2)
760 ReturnVal = ReturnVal2;
761 }
762 else {
763 ReturnVal = this->Multiply(false, X, *T, 1.0, 0.0);
764 int ReturnVal2 = this->Multiply(true, *T, Y, Alpha, Beta);
765 if(ReturnVal < ReturnVal2)
766 ReturnVal = ReturnVal2;
767 }
768 }
769 }
770 if (Xcopy!=0) {
771 delete Xcopy;
772 EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of x
773 return(1);
774 }
775 return ReturnVal;
776}
777
779 const Epetra_Vector& x,
780 Epetra_Vector& y,
781 const Epetra_Vector& w,
782 Epetra_Vector& z,
783 double Alpha,
784 double Beta,
785 double Omega,
786 double Zeta) const {
787
788 int ReturnVal;
789
790 if(!Filled())
791 EPETRA_CHK_ERR(-1); // Matrix must be filled.
792
793 double* xp = (double*) x.Values();
794 double* xp2 = (double*) x.Values();
795 double* wp = (double*) w.Values();
796 double* yp = (double*) y.Values();
797// double* yp2 = (double*) y.Values();
798 double* zp = (double*) z.Values();
799 Epetra_MultiVector* yp2 = 0;
800 Epetra_Vector* xcopy = 0;
801 if (&x==&y && Importer()==0 && Exporter()==0) {
802 xcopy = new Epetra_Vector(x);
803 xp = (double *) xcopy->Values();
804 }
805 Epetra_Vector* wcopy = 0;
806 if (&w==&z && Importer()==0 && Exporter()==0) {
807 wcopy = new Epetra_Vector(w);
808 wp = (double *) wcopy->Values();
809 }
810 UpdateImportVector(1); // Refresh import and output vectors if needed
812
813 if(TransA) {
814
815 if(Importer() != 0) {
817 xp = (double*) ImportVector_->Values();
818 xp2 = new double[ImportVector_->MyLength()];
819 for(int i = 0; i < ImportVector_->MyLength(); i++)
820 xp2[i] = xp[i];
822 zp = (double*) ImportVector_->Values();
823 }
824 if(Exporter() != 0) {
826 yp = (double*) yp2->Values();
827
828 //for(int i = 0; i < ExportVector_->MyLength(); i++)
829 //yp2[i] = yp[i];
830 wp = (double*) ExportVector_->Values();
831 }
832
833 oski_vecview_t oskiX=0;
834 oski_vecview_t oskiY=0;
835 oski_vecview_t oskiW=0;
836 oski_vecview_t oskiZ=0;
837 if(Importer() != 0) {
838 oskiX = oski_CreateVecView(xp2,ImportVector_->MyLength(),1);
839 oskiZ = oski_CreateVecView(zp,ImportVector_->MyLength(),1);
840 }
841 else {
842 oskiX = oski_CreateVecView(xp2,x.MyLength(),1);
843 oskiZ = oski_CreateVecView(zp,z.MyLength(),1);
844 }
845 if(Exporter() != 0) {
846 oskiY = oski_CreateVecView(yp,ExportVector_->MyLength(),1);
847 oskiW = oski_CreateVecView(wp,ExportVector_->MyLength(),1);
848 }
849 else {
850 oskiY = oski_CreateVecView(yp,y.MyLength(),1);
851 oskiW = oski_CreateVecView(wp,w.MyLength(),1);
852 }
853
854 ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, oskiX, Beta, oskiY, OP_TRANS, Omega, oskiW, Zeta, oskiZ);
855 if(Exporter() != 0) {
856 y.PutScalar(0.0); // Make sure target is zero
857 EPETRA_CHK_ERR(y.Export(*yp2, *Exporter(), Add)); // Fill y with Values from export vector
858 }
859 if(Importer() != 0) {
860 z.PutScalar(0.0); // Make sure target is zero
861 EPETRA_CHK_ERR(z.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
862 }
863 // Handle case of rangemap being a local replicated map
864 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
865 if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(z.Reduce());
866
868 }
869 // ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, (*xCast).Oski_View(), Beta, (*yCast).Oski_View(), OP_TRANS, Omega, (*wCast).Oski_View(), Zeta, (*zCast).Oski_View());
870 else {
871
872 if(Importer() != 0) {
874 xp = (double*) ImportVector_->Values();
875 xp2 = new double[ImportVector_->MyLength()];
876 for(int i = 0; i < ImportVector_->MyLength(); i++)
877 xp2[i] = xp[i];
879 wp = (double*) ImportVector_->Values();
880 }
881 if(Exporter() != 0) {
883 yp = (double*) yp2->Values();
884
885 //for(int i = 0; i < ExportVector_->MyLength(); i++)
886 //yp2[i] = yp[i];
887 zp = (double*) ExportVector_->Values();
888 }
889
890 oski_vecview_t oskiX=0;
891 oski_vecview_t oskiY=0;
892 oski_vecview_t oskiW=0;
893 oski_vecview_t oskiZ=0;
894 if(Importer() != 0) {
895 oskiX = oski_CreateVecView(xp2,ImportVector_->MyLength(),1);
896 oskiW = oski_CreateVecView(wp,ImportVector_->MyLength(),1);
897 }
898 else {
899 oskiX = oski_CreateVecView(xp2,x.MyLength(),1);
900 oskiW = oski_CreateVecView(wp,w.MyLength(),1);
901 }
902 if(Exporter() != 0) {
903 oskiY = oski_CreateVecView(yp,ExportVector_->MyLength(),1);
904 oskiZ = oski_CreateVecView(zp,ExportVector_->MyLength(),1);
905 }
906 else {
907 oskiY = oski_CreateVecView(yp,y.MyLength(),1);
908 oskiZ = oski_CreateVecView(zp,z.MyLength(),1);
909 }
910
911 ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, oskiX, Beta, oskiY, OP_NORMAL, Omega, oskiW, Zeta, oskiZ);
912 if(Exporter() != 0) {
913 y.PutScalar(0.0); // Make sure target is zero
914 EPETRA_CHK_ERR(y.Export(*yp2, *Exporter(), Add)); // Fill y with Values from export vector
915 z.PutScalar(0.0); // Make sure target is zero
916 EPETRA_CHK_ERR(z.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector*/
917 }
918 // Handle case of rangemap being a local replicated map
919 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
920 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(z.Reduce());
921
923
924 }
925 if(ReturnVal)
926 std::cerr << "OskiVector multiply error\n";
927 return ReturnVal;
928}
929
931 const Epetra_MultiVector& X,
933 const Epetra_MultiVector& W,
935 double Alpha,
936 double Beta,
937 double Omega,
938 double Zeta) const {
939 int ReturnVal;
940 if(!Filled())
941 EPETRA_CHK_ERR(-1); // Matrix must be filled.
942 int NumVectors = X.NumVectors();
943 if (NumVectors!=Y.NumVectors()) {
944 EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
945 }
946
947
948 double** Xp = (double**) X.Pointers();
949 double** Xp2 = (double**) X.Pointers();
950 double** Wp = (double**) W.Pointers();
951 double** Yp = (double**) Y.Pointers();
952 double** Zp = (double**) Z.Pointers();
953 int LDX = X.ConstantStride() ? X.Stride() : 0;
954 int LDY = Y.ConstantStride() ? Y.Stride() : 0;
955 int LDW = W.ConstantStride() ? W.Stride() : 0;
956 int LDZ = Z.ConstantStride() ? Z.Stride() : 0;
957
958 Epetra_MultiVector* Yp2 = 0;
959 Epetra_MultiVector* X2 = 0;
960 Epetra_MultiVector* Xcopy = 0;
961 if (&X==&Y && Importer()==0 && Exporter()==0) {
962 Xcopy = new Epetra_MultiVector(X);
963 Xp = (double **) Xcopy->Pointers();
964 LDX = Xcopy->ConstantStride() ? Xcopy->Stride() : 0;
965 }
966 Epetra_MultiVector* Wcopy = 0;
967 if (&W==&Z && Importer()==0 && Exporter()==0) {
968 Wcopy = new Epetra_MultiVector(W);
969 Wp = (double **) Wcopy->Values();
970 LDW = Wcopy->ConstantStride() ? Wcopy->Stride() : 0;
971 }
972 UpdateImportVector(NumVectors); // Refresh import and output vectors if needed
973 UpdateExportVector(NumVectors);
974
975 if(TransA) {
976
977 if(Importer() != 0) {
979 Xp = (double**) ImportVector_->Pointers();
981 X2 = new Epetra_MultiVector(X);
982 double** Xp2 = (double**) X2->Pointers();
983 Xp2 = (double **) X2->Pointers();
985 Zp = (double**) ImportVector_->Pointers();
987 }
988 if(Exporter() != 0) {
990 Yp = (double**) Yp2->Pointers();
991 Wp = (double**) ExportVector_->Pointers();
992 }
993
994 oski_vecview_t oskiX=0;
995 oski_vecview_t oskiY=0;
996 oski_vecview_t oskiW=0;
997 oski_vecview_t oskiZ=0;
998 if(Importer() != 0) {
999 oskiX = oski_CreateMultiVecView(*Xp2,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
1000 oskiZ = oski_CreateMultiVecView(*Zp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDZ);
1001 }
1002 else {
1003 oskiX = oski_CreateMultiVecView(*Xp2,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
1004 oskiZ = oski_CreateMultiVecView(*Zp,Z.MyLength(),NumVectors,LAYOUT_COLMAJ,LDZ);
1005 }
1006 if(Exporter() != 0) {
1007 oskiY = oski_CreateMultiVecView(*Yp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
1008 oskiW = oski_CreateMultiVecView(*Wp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDW);
1009 }
1010 else {
1011 oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
1012 oskiW = oski_CreateMultiVecView(*Wp,W.MyLength(),NumVectors,LAYOUT_COLMAJ,LDW);
1013 }
1014
1015 ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, oskiX, Beta, oskiY, OP_TRANS, Omega, oskiW, Zeta, oskiZ);
1016 if(Exporter() != 0) {
1017 Y.PutScalar(0.0); // Make sure target is zero
1018 EPETRA_CHK_ERR(Y.Export(*Yp2, *Exporter(), Add)); // Fill y with Values from export vector
1019 }
1020 if(Importer() != 0) {
1021 Z.PutScalar(0.0); // Make sure target is zero
1022 EPETRA_CHK_ERR(Z.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
1023 }
1024 // Handle case of rangemap being a local replicated map
1025 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
1026 if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Z.Reduce());
1027
1029 }
1030 // ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, (*xCast).Oski_View(), Beta, (*yCast).Oski_View(), OP_TRANS, Omega, (*wCast).Oski_View(), Zeta, (*zCast).Oski_View());
1031 else {
1032
1033 if(Importer() != 0) {
1035 Xp = (double**) ImportVector_->Pointers();
1037 X2 = new Epetra_MultiVector(X);
1038 Xp2 = (double**) X2->Pointers();
1040 Wp = (double**) ImportVector_->Pointers();
1042 }
1043 if(Exporter() != 0) {
1045 Yp = (double**) Yp2->Pointers();
1046 Zp = (double**) ExportVector_->Pointers();
1047 }
1048
1049 oski_vecview_t oskiX=0;
1050 oski_vecview_t oskiY=0;
1051 oski_vecview_t oskiW=0;
1052 oski_vecview_t oskiZ=0;
1053 if(Importer() != 0) {
1054 oskiX = oski_CreateMultiVecView(*Xp2,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
1055 oskiW = oski_CreateMultiVecView(*Wp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDW);
1056 }
1057 else {
1058 oskiX = oski_CreateMultiVecView(*Xp2,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
1059 oskiW = oski_CreateMultiVecView(*Wp,W.MyLength(),NumVectors,LAYOUT_COLMAJ,LDW);
1060 }
1061 if(Exporter() != 0) {
1062 oskiY = oski_CreateMultiVecView(*Yp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
1063 oskiZ = oski_CreateMultiVecView(*Zp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDZ);
1064 }
1065 else {
1066 oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
1067 oskiZ = oski_CreateMultiVecView(*Zp,Z.MyLength(),NumVectors,LAYOUT_COLMAJ,LDZ);
1068 }
1069
1070 ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, oskiX, Beta, oskiY, OP_NORMAL, Omega, oskiW, Zeta, oskiZ);
1071 if(Exporter() != 0) {
1072 Y.PutScalar(0.0); // Make sure target is zero
1073 EPETRA_CHK_ERR(Y.Export(*Yp2, *Exporter(), Add)); // Fill y with Values from export vector
1074 Z.PutScalar(0.0); // Make sure target is zero
1075 EPETRA_CHK_ERR(Z.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector*/
1076 }
1077 // Handle case of rangemap being a local replicated map
1078 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
1079 if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Z.Reduce());
1080
1082
1083 }
1084 return ReturnVal;
1085}
1086
1087int Epetra_OskiMatrix::MatPowMultiply(bool TransA,
1088 const Epetra_Vector& x,
1089 Epetra_Vector& y,
1091 int Power,
1092 double Alpha,
1093 double Beta) const {
1094 //The code has not been tested. It should work in serial but not in parallel.
1095 std::cerr << "MatPowMultiply not implemented in oski-1.01h release.\n";
1096 return -1;
1097
1098 int ReturnVal;
1099
1100 if(!Filled())
1101 EPETRA_CHK_ERR(-1); // Matrix must be filled.
1102
1103 double* xp = (double*) x.Values();
1104 double* yp = (double*) y.Values();
1105 double** Tp = (double**) T.Pointers();
1106
1107 Epetra_MultiVector *Tptr;
1108
1109 int LDT = T.ConstantStride() ? T.Stride() : 0;
1110
1111
1112 if(this->Comm().NumProc() == 1) {
1113 oski_vecview_t oskiX=0;
1114 oski_vecview_t oskiY=0;
1115 oski_vecview_t oskiT=0;
1116 if (&T != NULL)
1117 oskiT = oski_CreateMultiVecView(*Tp,T.MyLength(),Power,LAYOUT_COLMAJ,LDT);
1118 oskiX = oski_CreateVecView(xp,x.MyLength(),1);
1119 oskiY = oski_CreateVecView(yp,y.MyLength(),1);
1120
1121 if(TransA)
1122 ReturnVal = oski_MatPowMult(A_tunable_, OP_TRANS, Power, Alpha, oskiX, Beta, oskiY, oskiT);
1123 else
1124 ReturnVal = oski_MatPowMult(A_tunable_, OP_NORMAL, Power, Alpha, oskiX, Beta, oskiY, oskiT);
1125 if(ReturnVal)
1126 std::cerr << "OskiVector matpow multiply error\n";
1127 }
1128 else {
1129 if(&T == NULL)
1130 Tptr = new Epetra_MultiVector(x.Map(), Power);
1131 else
1132 Tptr = &T;
1133 if(TransA) {
1134 ReturnVal = this->Multiply(true, x, Tptr[1], 1.0, 0.0);
1135 for(int i = 1; i < Power-1; i++)
1136 ReturnVal = this->Multiply(true, Tptr[i], Tptr[i+1], 1.0, 0.0);
1137 ReturnVal = this->Multiply(true, Tptr[Power-2], y, Alpha, Beta);
1138 }
1139 else {
1140 ReturnVal = this->Multiply(false, x, Tptr[1], 1.0, 0.0);
1141 for(int i = 1; i < Power-1; i++)
1142 ReturnVal = this->Multiply(false, Tptr[i], Tptr[i+1], 1.0, 0.0);
1143 ReturnVal = this->Multiply(false, Tptr[Power-2], y, Alpha, Beta);
1144 }
1145 if(ReturnVal)
1146 std::cerr << "OskiVector matpow multiply error\n";
1147 if(&T == NULL)
1148 delete(Tptr);
1149 }
1150 UpdateFlops(2 * Power * NumGlobalNonzeros64());
1151 return ReturnVal;
1152}
1153
1154int Epetra_OskiMatrix::MatPowMultiply(bool TransA,
1155 const Epetra_Vector& x,
1156 Epetra_Vector& y,
1157 int Power,
1158 double Alpha,
1159 double Beta) const {
1160
1161 //The code has not been tested. It should work in serial but not in parallel.
1162 std::cerr << "MatPowMultiply not implemented in oski-1.01h release.\n";
1163 return -1;
1164 std::cerr << "mult\n";
1165 Epetra_OskiVector* xCast = NULL;
1166 Epetra_OskiVector* yCast = NULL;
1167 bool xCreate = false;
1168 bool yCreate = false;
1169 int ReturnVal;
1170 xCast = dynamic_cast<Epetra_OskiVector*>(const_cast <Epetra_Vector*>(&x));
1171 yCast = dynamic_cast<Epetra_OskiVector*>(&y);
1172 if (xCast == NULL) {
1173 xCast = new Epetra_OskiVector(x);
1174 xCreate = true;
1175 }
1176 if (yCast == NULL) {
1177 yCast = new Epetra_OskiVector(y);
1178 yCreate = true;
1179 }
1180 std::cerr << "mult\n";
1181 if(TransA)
1182 ReturnVal = oski_MatPowMult(A_tunable_, OP_TRANS, Power, Alpha, (*xCast).Oski_View(), Beta, (*yCast).Oski_View(), NULL);
1183 else
1184 ReturnVal = oski_MatPowMult(A_tunable_, OP_NORMAL, Power, Alpha, (*xCast).Oski_View(), Beta, (*yCast).Oski_View(), NULL);
1185 std::cerr << "done\n";
1186 if(ReturnVal)
1187 std::cerr << "OskiVector matpow multiply error\n";
1188 if(xCreate)
1189 delete xCast;
1190 if(yCreate)
1191 delete yCast;
1192 return ReturnVal;
1193}
1194
1195int Epetra_OskiMatrix::SetHint(const Teuchos::ParameterList& List) {
1196 int* ArgArray = NULL;
1197 int Diags;
1198 int Blocks;
1199 char Number[10];
1200 char Row[20];
1201 char Col[20];
1202 char Diag[20];
1203 int ReturnVal = 0;
1204 if(List.isParameter("randompattern"))
1205 if(Teuchos::getParameter<bool>(List, "randompattern"))
1206 if(ReturnVal = oski_SetHint(A_tunable_, HINT_RANDOM_PATTERN))
1207 std::cerr << "Error setting hint random pattern.\n";
1208 if(List.isParameter("correlatedpattern"))
1209 if(Teuchos::getParameter<bool>(List, "correlatedpattern"))
1210 if(ReturnVal = oski_SetHint(A_tunable_, HINT_CORREL_PATTERN))
1211 std::cerr << "Error setting hint correlated pattern.\n";
1212 if(List.isParameter("symmetricpattern"))
1213 if(Teuchos::getParameter<bool>(List, "symmetricpattern"))
1214 if(ReturnVal = oski_SetHint(A_tunable_, HINT_SYMM_PATTERN))
1215 std::cerr << "Error setting hint symmetric pattern.\n";
1216 if(List.isParameter("nonsymmetricpattern"))
1217 if(Teuchos::getParameter<bool>(List, "nonsymmetricpattern"))
1218 if(ReturnVal = oski_SetHint(A_tunable_, HINT_NONSYMM_PATTERN))
1219 std::cerr << "Error setting hint nonsymmetric pattern.\n";
1220 if(List.isParameter("alignedblocks"))
1221 if(Teuchos::getParameter<bool>(List, "alignedblocks"))
1222 {
1223 if(ReturnVal = oski_SetHint(A_tunable_, HINT_ALIGNED_BLOCKS))
1224 std::cerr << "Error setting hint aligned blocks.\n";
1225 }
1226 if(List.isParameter("unalignedblocks"))
1227 if(Teuchos::getParameter<bool>(List, "unalignedblocks"))
1228 if(ReturnVal = oski_SetHint(A_tunable_, HINT_UNALIGNED_BLOCKS))
1229 std::cerr << "Error setting hint unaligned blocks.\n";
1230 if(List.isParameter("nodiags"))
1231 if(Teuchos::getParameter<bool>(List, "nodiags"))
1232 if(ReturnVal = oski_SetHint(A_tunable_, HINT_NO_DIAGS))
1233 std::cerr << "Error setting hint no diags.\n";
1234 if(List.isParameter("noblocks"))
1235 if(Teuchos::getParameter<bool>(List, "noblocks"))
1236 if(ReturnVal = oski_SetHint(A_tunable_, HINT_NO_BLOCKS))
1237 std::cerr << "Error setting hint no blocks.\n";
1238 if(List.isParameter("singleblocksize"))
1239 if(Teuchos::getParameter<bool>(List, "singleblocksize")) {
1240 ArgArray = new int[2];
1241 if(List.isParameter("row"))
1242 {
1243 ArgArray[0] = Teuchos::getParameter<int>(List, "row");
1244 if(List.isParameter("col"))
1245 ArgArray[1] = Teuchos::getParameter<int>(List, "col");
1246 if(ReturnVal = oski_SetHint(A_tunable_, HINT_SINGLE_BLOCKSIZE, ArgArray[0], ArgArray[1]))
1247 std::cerr << "Error setting hint single block size.\n";
1248 }
1249 else
1250 if(ReturnVal = oski_SetHint(A_tunable_, HINT_SINGLE_BLOCKSIZE))
1251 std::cerr << "Error setting hint single block size.\n";
1252 delete [] ArgArray;
1253 ArgArray = NULL;
1254 }
1255 if(List.isParameter("multipleblocksize"))
1256 if(Teuchos::getParameter<bool>(List, "multipleblocksize"))
1257 if(List.isParameter("blocks")) {
1258 Blocks = Teuchos::getParameter<int>(List, "blocks");
1259 ArgArray = new int[2*Blocks+1];
1260 ArgArray[0] = Blocks;
1261 for(int i = 0; i < Blocks; i++) {
1262 sprintf(Number, "%d", i+1);
1263 strcpy(Row, "row");
1264 strcpy(Col, "col");
1265 strcat(Row, Number);
1266 strcat(Col, Number);
1267 if(List.isParameter(Row))
1268 ArgArray[i*2 + 1] = Teuchos::getParameter<int>(List, Row);
1269 if(List.isParameter(Col))
1270 ArgArray[i*2 + 2] = Teuchos::getParameter<int>(List, Col);
1271 }
1272 switch(Blocks) {
1273 case 1 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2]))
1274 std::cerr << "Error setting hint multiple blocks.\n";
1275 break;
1276 case 2 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4]))
1277 std::cerr << "Error setting hint multiple blocks.\n";
1278 break;
1279 case 3 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4], ArgArray[5], ArgArray[6]))
1280 std::cerr << "Error setting hint multiple blocks.\n";
1281 break;
1282 case 4 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4], ArgArray[5], ArgArray[6], ArgArray[7], ArgArray[8]))
1283 std::cerr << "Error setting hint multiple blocks.\n";
1284 break;
1285 case 5 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4], ArgArray[5], ArgArray[6], ArgArray[7], ArgArray[8], ArgArray[9], ArgArray[10]))
1286 std::cerr << "Error setting hint multiple blocks.\n";
1287 break;
1288 default : if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES))
1289 std::cerr << "Error setting hint multiple blocks.\n";
1290 break;
1291 }
1292 delete [] ArgArray;
1293 ArgArray = NULL;
1294 }
1295 else
1296 if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES))
1297 std::cerr << "Error setting hint multiple blocks.\n";
1298 if(List.isParameter("diags"))
1299 if(Teuchos::getParameter<bool>(List, "diags"))
1300 if(List.isParameter("numdiags")) {
1301 Diags = Teuchos::getParameter<int>(List, "numdiags");
1302 ArgArray = new int[Diags+1];
1303 ArgArray[0] = Diags;
1304 for(int i = 0; i < Diags; i++) {
1305 sprintf(Number, "%d", i + 1);
1306 strcpy(Diag, "diag");
1307 strcat(Diag, Number);
1308 if(List.isParameter(Diag))
1309 ArgArray[i + 1] = Teuchos::getParameter<int>(List, Diag);
1310 }
1311 switch(Diags) {
1312 case 1 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1]))
1313 std::cerr << "Error setting hint diags\n";
1314 break;
1315 case 2 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1], ArgArray[2]))
1316 std::cerr << "Error setting hint diags\n";
1317 break;
1318 case 3 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3]))
1319 std::cerr << "Error setting hint diags\n";
1320 break;
1321 case 4 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4]))
1322 std::cerr << "Error setting hint diags\n";
1323 break;
1324 case 5 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4], ArgArray[5]))
1325 std::cerr << "Error setting hint diags\n";
1326 break;
1327 default : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0]))
1328 std::cerr << "Error setting hint diags\n";
1329 break;
1330 }
1331 delete [] ArgArray;
1332 }
1333 else
1334 {
1335 if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS))
1336 std::cerr << "Error setting hint digs.\n";
1337 }
1338 return ReturnVal;
1339}
1340
1342 double Alpha,
1343 const Epetra_OskiMultiVector& InVec,
1344 double Beta,
1345 const Epetra_OskiMultiVector& OutVec,
1346 int NumCalls,
1347 const Teuchos::ParameterList& List) {
1348 int ReturnVal;
1349 oski_vecview_t InView = NULL;
1350 oski_vecview_t OutView = NULL;
1351 InView = InVec.Oski_View();
1352 OutView = OutVec.Oski_View();
1353 if(List.isParameter("tune"))
1354 if(Teuchos::getParameter<bool>(List, "tune"))
1355 NumCalls = ALWAYS_TUNE;
1356 if(List.isParameter("tuneaggressive"))
1357 if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
1358 NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
1359 if(List.isParameter("symminvec"))
1360 if(Teuchos::getParameter<bool>(List, "symminvec"))
1361 InView = SYMBOLIC_VEC;
1362 if(List.isParameter("symminmultivec"))
1363 if(Teuchos::getParameter<bool>(List, "symminmultivec"))
1364 InView = SYMBOLIC_MULTIVEC;
1365 if(List.isParameter("symmoutvec"))
1366 if(Teuchos::getParameter<bool>(List, "symmoutvec"))
1367 OutView = SYMBOLIC_VEC;
1368 if(List.isParameter("symmoutmultivec"))
1369 if(Teuchos::getParameter<bool>(List, "symmoutmultivec"))
1370 OutView = SYMBOLIC_MULTIVEC;
1371 if(this->Comm().NumProc() > 1) {
1372 if(InVec.NumVectors() == 1)
1373 InView = SYMBOLIC_VEC;
1374 else
1375 InView = SYMBOLIC_MULTIVEC;
1376 if(OutVec.NumVectors() == 1)
1377 OutView = SYMBOLIC_VEC;
1378 else
1379 OutView = SYMBOLIC_MULTIVEC;
1380 }
1381 if(TransA)
1382 ReturnVal = oski_SetHintMatMult(A_tunable_, OP_TRANS, Alpha, InView, Beta, OutView, NumCalls);
1383 else
1384 ReturnVal = oski_SetHintMatMult(A_tunable_, OP_NORMAL, Alpha, InView, Beta, OutView, NumCalls);
1385 if(ReturnVal)
1386 std::cerr << "Set hint multivector multiply error\n";
1387 return ReturnVal;
1388}
1389
1390int Epetra_OskiMatrix::SetHintSolve(bool TransA,
1391 double Alpha,
1392 const Epetra_OskiMultiVector& Vector,
1393 int NumCalls,
1394 const Teuchos::ParameterList& List) {
1395 int ReturnVal;
1396 oski_vecview_t VecView = NULL;
1397 VecView = Vector.Oski_View();
1398 if(List.isParameter("tune"))
1399 if(Teuchos::getParameter<bool>(List, "tune"))
1400 NumCalls = ALWAYS_TUNE;
1401 if(List.isParameter("tuneaggressive"))
1402 if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
1403 NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
1404 if(List.isParameter("symmvec"))
1405 if(Teuchos::getParameter<bool>(List, "symmvec"))
1406 VecView = SYMBOLIC_VEC;
1407 if(List.isParameter("symmmultivec"))
1408 if(Teuchos::getParameter<bool>(List, "symmmultivec"))
1409 VecView = SYMBOLIC_MULTIVEC;
1410 if(this->Comm().NumProc() > 1) {
1411 if(Vector.NumVectors() == 1)
1412 VecView = SYMBOLIC_VEC;
1413 else
1414 VecView = SYMBOLIC_MULTIVEC;
1415 }
1416 if(TransA)
1417 ReturnVal = oski_SetHintMatTrisolve(A_tunable_, OP_TRANS, Alpha, VecView, NumCalls);
1418 else
1419 ReturnVal = oski_SetHintMatTrisolve(A_tunable_, OP_NORMAL, Alpha, VecView, NumCalls);
1420 if(ReturnVal)
1421 std::cerr << "Set hint solve error\n";
1422 return ReturnVal;
1423}
1424
1426 double Alpha,
1427 const Epetra_OskiMultiVector& InVec,
1428 double Beta,
1429 const Epetra_OskiMultiVector& OutVec,
1430 const Epetra_OskiMultiVector& Intermediate,
1431 int NumCalls,
1432 const Teuchos::ParameterList& List) {
1433 int ReturnVal;
1434 oski_vecview_t InView = NULL;
1435 oski_vecview_t OutView = NULL;
1436 oski_vecview_t IntermediateView = NULL;
1437 InView = InVec.Oski_View();
1438 OutView = OutVec.Oski_View();
1439 if(&Intermediate != NULL)
1440 IntermediateView = Intermediate.Oski_View();
1441 if(List.isParameter("tune"))
1442 if(Teuchos::getParameter<bool>(List, "tune"))
1443 NumCalls = ALWAYS_TUNE;
1444 if(List.isParameter("tuneaggressive"))
1445 if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
1446 NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
1447 if(List.isParameter("symminvec"))
1448 if(Teuchos::getParameter<bool>(List, "symminvec"))
1449 InView = SYMBOLIC_VEC;
1450 if(List.isParameter("symminmultivec"))
1451 if(Teuchos::getParameter<bool>(List, "symminmultivec"))
1452 InView = SYMBOLIC_MULTIVEC;
1453 if(List.isParameter("symmoutvec"))
1454 if(Teuchos::getParameter<bool>(List, "symmoutvec"))
1455 OutView = SYMBOLIC_VEC;
1456 if(List.isParameter("symmoutmultivec"))
1457 if(Teuchos::getParameter<bool>(List, "symmoutmultivec"))
1458 OutView = SYMBOLIC_MULTIVEC;
1459 if(List.isParameter("symmintervec"))
1460 if(Teuchos::getParameter<bool>(List, "symmintervec"))
1461 IntermediateView = SYMBOLIC_VEC;
1462 if(List.isParameter("symmintermultivec"))
1463 if(Teuchos::getParameter<bool>(List, "symmintermultivec"))
1464 IntermediateView = SYMBOLIC_MULTIVEC;
1465 if(List.isParameter("invalidinter"))
1466 if(Teuchos::getParameter<bool>(List, "invalidinter"))
1467 IntermediateView = NULL;
1468 if(this->Comm().NumProc() > 1) {
1469 if(InVec.NumVectors() == 1)
1470 InView = SYMBOLIC_VEC;
1471 else
1472 InView = SYMBOLIC_MULTIVEC;
1473 if(OutVec.NumVectors() == 1)
1474 OutView = SYMBOLIC_VEC;
1475 else
1476 OutView = SYMBOLIC_MULTIVEC;
1477 if(Intermediate.NumVectors() == 1)
1478 IntermediateView = SYMBOLIC_VEC;
1479 else
1480 IntermediateView = SYMBOLIC_MULTIVEC;
1481 }
1482 if(ATA)
1483 ReturnVal = oski_SetHintMatTransMatMult(A_tunable_, OP_AT_A, Alpha, InView, Beta, OutView, IntermediateView, NumCalls);
1484 else
1485 ReturnVal = oski_SetHintMatTransMatMult(A_tunable_, OP_A_AT, Alpha, InView, Beta, OutView, IntermediateView, NumCalls);
1486 if(ReturnVal)
1487 std::cerr << "Set hint mattransmat multiply error\n";
1488 return ReturnVal;
1489}
1490
1492 double Alpha,
1493 const Epetra_OskiMultiVector& InVec,
1494 double Beta,
1495 const Epetra_OskiMultiVector& OutVec,
1496 double Omega,
1497 const Epetra_OskiMultiVector& InVec2,
1498 double Zeta,
1499 const Epetra_OskiMultiVector& OutVec2,
1500 int NumCalls,
1501 const Teuchos::ParameterList& List) {
1502 int ReturnVal;
1503 oski_vecview_t InView = NULL;
1504 oski_vecview_t OutView = NULL;
1505 oski_vecview_t InView2 = NULL;
1506 oski_vecview_t OutView2 = NULL;
1507 InView = InVec.Oski_View();
1508 OutView = OutVec.Oski_View();
1509 InView2 = InVec2.Oski_View();
1510 OutView2 = OutVec2.Oski_View();
1511 if(List.isParameter("tune"))
1512 if(Teuchos::getParameter<bool>(List, "tune"))
1513 NumCalls = ALWAYS_TUNE;
1514 if(List.isParameter("tuneaggressive"))
1515 if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
1516 NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
1517 if(List.isParameter("symminvec"))
1518 if(Teuchos::getParameter<bool>(List, "symminvec"))
1519 InView = SYMBOLIC_VEC;
1520 if(List.isParameter("symminmultivec"))
1521 if(Teuchos::getParameter<bool>(List, "symminmultivec"))
1522 InView = SYMBOLIC_MULTIVEC;
1523 if(List.isParameter("symmoutvec"))
1524 if(Teuchos::getParameter<bool>(List, "symmoutvec"))
1525 OutView = SYMBOLIC_VEC;
1526 if(List.isParameter("symmoutmultivec"))
1527 if(Teuchos::getParameter<bool>(List, "symmoutmultivec"))
1528 OutView = SYMBOLIC_MULTIVEC;
1529 if(List.isParameter("symminvec2"))
1530 if(Teuchos::getParameter<bool>(List, "symminvec2"))
1531 InView2 = SYMBOLIC_VEC;
1532 if(List.isParameter("symminmultivec2"))
1533 if(Teuchos::getParameter<bool>(List, "symminmultivec2"))
1534 InView2 = SYMBOLIC_MULTIVEC;
1535 if(List.isParameter("symmoutvec2"))
1536 if(Teuchos::getParameter<bool>(List, "symmoutvec2"))
1537 OutView2 = SYMBOLIC_VEC;
1538 if(List.isParameter("symmoutmultivec2"))
1539 if(Teuchos::getParameter<bool>(List, "symmoutmultivec2"))
1540 OutView2 = SYMBOLIC_MULTIVEC;
1541 if(this->Comm().NumProc() > 1) {
1542 if(InVec.NumVectors() == 1)
1543 InView = SYMBOLIC_VEC;
1544 else
1545 InView = SYMBOLIC_MULTIVEC;
1546 if(OutVec.NumVectors() == 1)
1547 OutView = SYMBOLIC_VEC;
1548 else
1549 OutView = SYMBOLIC_MULTIVEC;
1550 if(InVec2.NumVectors() == 1)
1551 InView2 = SYMBOLIC_VEC;
1552 else
1553 InView2 = SYMBOLIC_MULTIVEC;
1554 if(OutVec2.NumVectors() == 1)
1555 OutView2 = SYMBOLIC_VEC;
1556 else
1557 OutView2 = SYMBOLIC_MULTIVEC;
1558 }
1559 if(TransA)
1560 ReturnVal = oski_SetHintMatMultAndMatTransMult(A_tunable_, Alpha, InView, Beta, OutView, OP_TRANS, Omega, InView2, Zeta, OutView2, NumCalls);
1561 else
1562 ReturnVal = oski_SetHintMatMultAndMatTransMult(A_tunable_, Alpha, InView, Beta, OutView, OP_NORMAL, Omega, InView2, Zeta, OutView2, NumCalls);
1563 if(ReturnVal)
1564 std::cerr << "Set hint multivector multiply and mattransmultiply error\n";
1565 return ReturnVal;
1566}
1567
1569 double Alpha,
1570 const Epetra_OskiMultiVector& InVec,
1571 double Beta,
1572 const Epetra_OskiMultiVector& OutVec,
1573 const Epetra_OskiMultiVector& Intermediate,
1574 int Power,
1575 int NumCalls,
1576 const Teuchos::ParameterList& List) {
1577 int ReturnVal;
1578 oski_vecview_t InView = NULL;
1579 oski_vecview_t OutView = NULL;
1580 oski_vecview_t IntermediateView = NULL;
1581 InView = InVec.Oski_View();
1582 OutView = OutVec.Oski_View();
1583 if(&Intermediate != NULL)
1584 IntermediateView = Intermediate.Oski_View();
1585 if(List.isParameter("tune"))
1586 if(Teuchos::getParameter<bool>(List, "tune"))
1587 NumCalls = ALWAYS_TUNE;
1588 if(List.isParameter("tuneaggressive"))
1589 if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
1590 NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
1591 if(List.isParameter("symminvec"))
1592 if(Teuchos::getParameter<bool>(List, "symminvec"))
1593 InView = SYMBOLIC_VEC;
1594 if(List.isParameter("symminmultivec"))
1595 if(Teuchos::getParameter<bool>(List, "symminmultivec"))
1596 InView = SYMBOLIC_MULTIVEC;
1597 if(List.isParameter("symmoutvec"))
1598 if(Teuchos::getParameter<bool>(List, "symmoutvec"))
1599 OutView = SYMBOLIC_VEC;
1600 if(List.isParameter("symmoutmultivec"))
1601 if(Teuchos::getParameter<bool>(List, "symmoutmultivec"))
1602 OutView = SYMBOLIC_MULTIVEC;
1603 if(List.isParameter("symmintervec"))
1604 if(Teuchos::getParameter<bool>(List, "symmintermultivec"))
1605 IntermediateView = SYMBOLIC_MULTIVEC;
1606 if(List.isParameter("invalidinter"))
1607 if(Teuchos::getParameter<bool>(List, "invalidinter"))
1608 IntermediateView = NULL;
1609 if(this->Comm().NumProc() > 1) {
1610 if(InVec.NumVectors() == 1)
1611 InView = SYMBOLIC_VEC;
1612 else
1613 InView = SYMBOLIC_MULTIVEC;
1614 if(OutVec.NumVectors() == 1)
1615 OutView = SYMBOLIC_VEC;
1616 else
1617 OutView = SYMBOLIC_MULTIVEC;
1618 if(Intermediate.NumVectors() == 1)
1619 IntermediateView = SYMBOLIC_VEC;
1620 else
1621 IntermediateView = SYMBOLIC_MULTIVEC;
1622 }
1623 if(TransA)
1624 ReturnVal = oski_SetHintMatPowMult(A_tunable_, OP_TRANS, Power, Alpha, InView, Beta, OutView, IntermediateView, NumCalls);
1625 else
1626 ReturnVal = oski_SetHintMatPowMult(A_tunable_, OP_NORMAL, Power, Alpha, InView, Beta, OutView, IntermediateView, NumCalls);
1627 if(ReturnVal)
1628 std::cerr << "Set hint matpow multiply error\n";
1629 return ReturnVal;
1630}
1631
1632
1634 int ReturnVal;
1635 ReturnVal = oski_TuneMat(A_tunable_);
1637 Copy_Created_ = true;
1638 return ReturnVal;
1639}
1640
1642 return oski_IsMatPermuted(A_tunable_);
1643}
1644
1646 //might need a more efficient way to do this
1647 Epetra_OskiMatrix* Transformed = NULL;
1648 Epetra_OskiMatrix Temp(*this); //should be some lightweight copy
1649 Transformed = &Temp;
1650 Transformed->A_tunable_ = const_cast <oski_matrix_t> (oski_ViewPermutedMat(A_tunable_));
1651 return *Transformed;
1652}
1653
1655 Epetra_OskiPermutation* RowPerm = NULL;
1656 RowPerm = new Epetra_OskiPermutation[1];
1657 (*RowPerm).ReplacePermutation(const_cast <oski_perm_t> (oski_ViewPermutedMatRowPerm(A_tunable_)));
1658 return *RowPerm;
1659}
1660
1662 Epetra_OskiPermutation* ColPerm;
1663 ColPerm = new Epetra_OskiPermutation[1];
1664 (*ColPerm).ReplacePermutation(const_cast <oski_perm_t> (oski_ViewPermutedMatColPerm(A_tunable_)));
1665 return *ColPerm;
1666}
1667
1669 char* ReturnVal = NULL;
1670 ReturnVal = oski_GetMatTransforms(A_tunable_);
1671 if(ReturnVal == NULL)
1672 std::cerr << "Error in GetMatrixTransforms\n";
1673 return ReturnVal;
1674}
1675
1676int Epetra_OskiMatrix::ApplyMatrixTransforms(const char* Transforms) {
1677 int ReturnVal;
1678 ReturnVal = oski_ApplyMatTransforms(A_tunable_, Transforms);
1679 if(ReturnVal)
1680 std::cerr << "Error in ApplyMatrixTransforms\n";
1681 return ReturnVal;
1682}
1683#endif
1684#endif
#define EPETRA_CHK_ERR(a)
virtual int NumProc() const =0
Returns total number of processes.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
const Epetra_Export * Exporter() const
Returns the Epetra_Export object that contains the export operations for distributed operations.
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
void UpdateImportVector(int NumVectors) const
void UpdateExportVector(int NumVectors) const
long long NumGlobalNonzeros64() const
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
const Epetra_Map & RangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
Epetra_MultiVector * ExportVector_
const Epetra_Map & DomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
const Epetra_Import * Importer() const
Returns the Epetra_Import object that contains the import operations for distributed operations.
double ** Values() const
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
Epetra_MultiVector * ImportVector_
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int NumVectors() const
Returns the number of vectors in the multi-vector.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true).
double * Values() const
Get pointer to MultiVector values.
double ** Pointers() const
Get pointer to individual vector pointers.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
Epetra_OskiMatrix: A class for constructing and using OSKI Matrices within Epetra....
int Solve(bool Upper, bool TransA, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
Performs a triangular solve of y = (this^TransA)^-1*x where this is a triangular matrix.
int IsMatrixTransformed() const
Returns 1 if the matrix has been reordered by tuning, and 0 if it has not been.
int ReplaceDiagonalValues(const Epetra_OskiVector &Diagonal)
Replaces diagonal values of the matrix with those in the user-provided vector.
Epetra_OskiMatrix(const Epetra_OskiMatrix &Source)
Copy constructor.
int SetHint(const Teuchos::ParameterList &List)
Stores the hints in List in the matrix structure.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Performs a matrix vector multiply of y = this^TransA*x.
int TuneMatrix()
Tunes the matrix multiply if its deemed profitable.
int MultiplyAndMatTransMultiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y, const Epetra_Vector &w, Epetra_Vector &z, double Alpha=1.0, double Beta=0.0, double Omega=1.0, double Zeta=0.0) const
Performs the two matrix vector multiplies of y = Alpha*this*x + Beta*y and z = Omega*this^TransA*w + ...
int SetHintMatTransMatMultiply(bool ATA, double Alpha, const Epetra_OskiMultiVector &InVec, double Beta, const Epetra_OskiMultiVector &OutVec, const Epetra_OskiMultiVector &Intermediate, int NumCalls, const Teuchos::ParameterList &List)
Workload hints for computing a two matrix-vector multiplies that are composed used by OskiTuneMat to ...
int ApplyMatrixTransforms(const char *Transforms)
Replaces the current data structure of the matrix with the one specified in Transforms.
const Epetra_OskiPermutation & ViewColumnPermutation() const
Returns a read only column/right permutation of the Matrix.
const Epetra_CrsMatrix * Epetra_View_
int SetHintMultiplyAndMatTransMultiply(bool TransA, double Alpha, const Epetra_OskiMultiVector &InVec, double Beta, const Epetra_OskiMultiVector &OutVec, double Omega, const Epetra_OskiMultiVector &InVec2, double Zeta, const Epetra_OskiMultiVector &OutVec2, int NumCalls, const Teuchos::ParameterList &List)
Workload hints for computing two matrix-vector multiplies used by OskiTuneMat to optimize the data st...
int ReplaceMyValues(int MyRow, int NumEntries, double *Values, int *Indices)
Replace current values with this list of entries for a given local row of the matrix....
int MatPowMultiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y, Epetra_MultiVector &T, int Power=2, double Alpha=1.0, double Beta=0.0) const
Performs a matrix vector multiply of y = Alpha*(this^TransA)^Power*x + Beta*y. This is not implemente...
int SetHintPowMultiply(bool TransA, double Alpha, const Epetra_OskiMultiVector &InVec, double Beta, const Epetra_OskiMultiVector &OutVec, const Epetra_OskiMultiVector &Intermediate, int Power, int NumCalls, const Teuchos::ParameterList &List)
Workload hints for computing a matrix-vector multiply performed Power times used by OskiTuneMat to op...
virtual ~Epetra_OskiMatrix()
Destructor
const Epetra_OskiPermutation & ViewRowPermutation() const
Returns a read only row/left permutation of the Matrix.
char * GetMatrixTransforms() const
Returns a string holding the transformations performed on the matrix when it was tuned.
const Epetra_OskiMatrix & ViewTransformedMat() const
Returns the transformed version of InMat if InMat has been transformed. If InMat has not been transfo...
int SetHintSolve(bool TransA, double Alpha, const Epetra_OskiMultiVector &Vector, int NumCalls, const Teuchos::ParameterList &List)
Workload hints for computing a triangular solve used by OskiTuneMat to optimize the data structure st...
int SumIntoMyValues(int MyRow, int NumEntries, double *Values, int *Indices)
Add this list of entries to existing values for a given local row of the matrix. WARNING: this could ...
int SetHintMultiply(bool TransA, double Alpha, const Epetra_OskiMultiVector &InVec, double Beta, const Epetra_OskiMultiVector &OutVec, int NumCalls, const Teuchos::ParameterList &List)
Workload hints for computing a matrix-vector multiply used by OskiTuneMat to optimize the data struct...
int MatTransMatMultiply(bool ATA, const Epetra_Vector &x, Epetra_Vector &y, Epetra_Vector *t, double Alpha=1.0, double Beta=0.0) const
Performs two matrix vector multiplies of y = Alpha*this^TransA*this*x + Beta*y or y = Alpha*this*this...
oski_matrix_t A_tunable_
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
Epetra_OskiMultiVector: A class for constructing and using dense Oski multi-vectors on a single proce...
oski_vecview_t Oski_View() const
Returns the Oski portion of the Multi-Vector.
Epetra_OskiPermutation: A class for storing the permutation performed on a Epetra_OskiMatrix.
void ReplacePermutation(const oski_perm_t &InPerm)
Stores a permutation in the data structure.
Epetra_OskiVector: A class for constructing and using dense OSKI vectors on a single processor or a s...
const Epetra_Vector * Epetra_View() const
Returns a view to the Epetra Object.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.