IFPACK Development
Loading...
Searching...
No Matches
Ifpack_SerialTriDiMatrix.cpp
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 "Ifpack_SerialTriDiMatrix.h"
44#include "Epetra_Util.h"
45
46//=============================================================================
49 Epetra_Object(-1, false),
50 N_(0),
51 A_Copied_(false),
52 CV_(Copy),
53 A_(0),
54 UseTranspose_(false)
55{
56 if (set_object_label) {
57 SetLabel("Epetra::SerialTriDiMatrix");
58 }
59}
60
61//=============================================================================
63 bool set_object_label)
65 Epetra_Object(-1, false),
66 N_(NumRowCol),
67 A_Copied_(false),
68 CV_(Copy),
69 A_(0),
70 UseTranspose_(false)
71{
72 if (set_object_label) {
73 SetLabel("Epetra::SerialTriDiMatrix");
74 }
75 if(NumRowCol < 0)
76 throw ReportError("NumRows = " + toString(NumRowCol) + ". Should be >= 0", -1);
77
78 int errorcode = Shape(NumRowCol);
79 if(errorcode != 0)
80 throw ReportError("Shape returned non-zero value", errorcode);
81}
82
83//=============================================================================
85 int NumRowCol,
86 bool set_object_label)
88 Epetra_Object(-1, false),
89 N_(NumRowCol),
90 A_Copied_(false),
91 CV_(CV_in),
92 A_(A_in),
93 UseTranspose_(false)
94{
95 if (set_object_label) {
96 SetLabel("Epetra::SerialTriDiMatrix");
97 }
98 if(A_in == 0)
99 throw ReportError("Null pointer passed as A parameter.", -3);
100 if(NumRowCol < 0)
101 throw ReportError("NumRowCol = " + toString(NumRowCol) + ". Should be >= 0", -1);
102
103 if (CV_in == Copy) {
104 const int newsize = (N_ == 1) ? 1 : 4*(N_-1);
105 if (newsize > 0) {
106 A_ = new double[newsize];
107 CopyMat(A_in, N_, A_, N_);
108 A_Copied_ = true;
109 DL_ = A_;
110 D_ = DL_+(N_-1);
111 DU_ = D_ + N_;
112 DU2_ = DU_ + (N_-1);
113 }
114 else {
115 A_ = 0;
116 }
117 }
118}
119//=============================================================================
121 : Epetra_CompObject(Source),
122 N_(Source.N_),
123 LDA_(Source.LDA_),
124 A_Copied_(false),
125 CV_(Source.CV_),
126 UseTranspose_(false)
127{
128 SetLabel(Source.Label());
129 if(CV_ == Copy) {
130 const int newsize = (N_ == 1)? 1 : 4*(N_-1);
131 if(newsize > 0) {
132 A_ = new double[newsize];
133 CopyMat(Source.A_, Source.N() , A_, N_);
134 A_Copied_ = true;
135 DL_ = A_;
136 D_ = DL_+(N_-1);
137 DU_ = D_ + N_;
138 DU2_ = DU_ + (N_-1);
139 }
140 else {
141 A_ = 0;
142 }
143 }
144}
145
146//=============================================================================
147int Ifpack_SerialTriDiMatrix::Reshape(int NumRows, int NumCols) {
148 if(NumRows < 0 || NumCols < 0)
149 return(-1);
150 if(NumRows != NumCols)
151 return(-1);
152
153 double* A_tmp = 0;
154 const int newsize = (N_ == 1)? 1 : 4*(N_-1);
155 if(newsize > 0) {
156 // Allocate space for new matrix
157 A_tmp = new double[newsize];
158 for (int k = 0; k < newsize; k++)
159 A_tmp[k] = 0.0; // Zero out values
160 if (A_ != 0)
161 CopyMat(A_, N_, A_tmp, NumRows); // Copy principal submatrix of A to new A
162
163 }
164 CleanupData(); // Get rid of anything that might be already allocated
165 N_ = NumCols;
166 if(newsize > 0) {
167 A_ = A_tmp; // Set pointer to new A
168 A_Copied_ = true;
169 }
170
171 DL_ = A_;
172 D_ = A_+(N_-1);
173 DU_ = D_+N_;
174 DU2_ = DU_+(N_-1);
175
176 return(0);
177}
178//=============================================================================
180 if(NumRowCol < 0 || NumRowCol < 0)
181 return(-1);
182
183 CleanupData(); // Get rid of anything that might be already allocated
184 N_ = NumRowCol;
185 LDA_ = N_;
186 const int newsize = (N_ == 1)? 1 : 4*(N_-1);
187 if(newsize > 0) {
188 A_ = new double[newsize];
189 for (int k = 0; k < newsize; k++)
190 A_[k] = 0.0; // Zero out values
191 A_Copied_ = true;
192 }
193 // set the pointers
194 DL_ = A_;
195 D_ = A_+(N_-1);
196 DU_ = D_+N_;
197 DU2_ = DU_+(N_-1);
198
199 return(0);
200}
201//=============================================================================
203{
204
205 CleanupData();
206
207}
208//=============================================================================
209void Ifpack_SerialTriDiMatrix::CleanupData()
210{
211 if (A_)
212 delete [] A_;
213 A_ = DL_ = D_ = DU_ = DU2_ = 0;
214 A_Copied_ = false;
215 N_ = 0;
216}
217//=============================================================================
219 if(this == &Source)
220 return(*this); // Special case of source same as target
221 if((CV_ == View) && (Source.CV_ == View) && (A_ == Source.A_))
222 return(*this); // Special case of both are views to same data.
223
224 if(std::strcmp(Label(), Source.Label()) != 0)
225 throw ReportError("operator= type mismatch (lhs = " + std::string(Label()) +
226 ", rhs = " + std::string(Source.Label()) + ").", -5);
227
228 if(Source.CV_ == View) {
229 if(CV_ == Copy) { // C->V only
230 CleanupData();
231 CV_ = View;
232 }
233 N_ = Source.N_;
234 A_ = Source.A_;
235 DL_ = Source.DL_;
236 D_ = Source.D_;
237 DU_ = Source.DU_;
238 DU2_ = Source.DU2_;
239 }
240 else {
241 if(CV_ == View) { // V->C
242 CV_ = Copy;
243 N_ = Source.N_;
244 const int newsize = 4*N_ - 4;
245 if(newsize > 0) {
246 A_ = new double[newsize];
247 DL_ = A_;
248 D_ = A_+(N_-1);
249 DU_ = D_+N_;
250 DU2_ = DU_+(N_-1);
251 A_Copied_ = true;
252 }
253 else {
254 A_ = 0;
255 DL_ = D_ = DU_ = DU2_ = 0;
256 A_Copied_ = false;
257 }
258 }
259 else { // C->C
260 if(Source.N_ == N_) { // we don't need to reallocate
261 N_ = Source.N_;
262 }
263 else { // we need to allocate more space (or less space)
264 CleanupData();
265 N_ = Source.N_;
266 const int newsize = (N_ == 1)? 1 : 4*(N_-1);
267 if(newsize > 0) {
268 A_ = new double[newsize];
269 DL_ = A_;
270 D_ = A_+(N_-1);
271 DU_ = D_+N_;
272 DU2_ = DU_+(N_-1);
273 A_Copied_ = true;
274 }
275 }
276 }
277 CopyMat(Source.A_, Source.N(), A_, N_); // V->C and C->C
278 }
279
280 return(*this);
281}
282
283
284//=============================================================================
286{
287 if (N_ != rhs.N_) return(false);
288
289 const double* A_tmp = A_;
290 const double* rhsA = rhs.A_;
291
292 const int size = (N_ == 1)? 1 : 4*(N_-1);
293
294 for(int j=0; j<size; ++j) {
295 if (std::abs(A_tmp[j] - rhsA[j]) > Epetra_MinDouble) {
296 return(false);
297 }
298 }
299
300 return(true);
301}
302
303//=============================================================================
305 if (N() != Source.N())
306 throw ReportError("Column dimension of source = " + toString(Source.N()) +
307 " is different than column dimension of target = " + toString(N()), -2);
308
309 CopyMat(Source.A(), Source.N(), A(), N(),true);
310 return(*this);
311}
312//=============================================================================
313void Ifpack_SerialTriDiMatrix::CopyMat(const double* Source,
314 int nrowcol,
315 double* Target,
316 int tN,
317 bool add)
318{
319 int lmax = EPETRA_MIN(nrowcol,tN);
320 if (add) {
321 // do this in 4 passes
322 for(int j=0; j<lmax; ++j) {
323 Target[(tN-1)+j] += Source[(nrowcol-1)+j]; //D
324 if(j<tN-1) {
325 Target[j] += Source[j]; // DL
326 Target[(tN-1)+tN + j] += Source[(nrowcol-1)+ nrowcol + j]; // DU
327 }
328 if(j<tN-2) Target[(tN-1)*2 + tN + j] += Source[ (nrowcol-1)*2 +nrowcol + j]; // DU2
329 }
330 }
331 else {
332 for(int j=0; j<lmax; ++j) {
333 Target[(tN-1)+j] = Source[(nrowcol-1)+j]; //D
334 if(j<tN-1) {
335 Target[j] = Source[j]; // DL
336 Target[(tN-1)+tN + j] = Source[(nrowcol-1)+ nrowcol + j]; // DU
337 }
338 if(j<tN-2) Target[(tN-1)*2 + tN + j] = Source[ (nrowcol-1)*2 +nrowcol + j]; // DU2
339 }
340 }
341 return;
342}
343//=============================================================================
345 int i;
346 double anorm = 0.0;
347 double * ptr = A_;
348 double sum=0.0;
349
350 const int size = (N_ == 1)? 1 : 4*(N_-1);
351
352 for (i=0; i<size; i++) sum += std::abs(*ptr++);
353
354 anorm = EPETRA_MAX(anorm, sum);
355 UpdateFlops((double)size );
356 return(anorm);
357}
358//=============================================================================
360 return NormOne();
361
362}
363//=============================================================================
365
366 int i;
367
368 double * ptr = A_;
369
370 const int size = (N_ == 1)? 1 : 4*(N_-1);
371
372 for (i=0; i<size ; i++) { *ptr = ScalarA * (*ptr); ptr++; }
373
374 UpdateFlops((double)N_*(double)N_);
375 return(0);
376
377}
378
379// //=========================================================================
380int Ifpack_SerialTriDiMatrix::Multiply (char /* TransA */, char /* TransB */, double /* ScalarAB */,
381 const Ifpack_SerialTriDiMatrix& /* A_in */,
382 const Ifpack_SerialTriDiMatrix& /* B */,
383 double /* ScalarThis */ ) {
384 throw ReportError("Ifpack_SerialTriDiMatrix::Multiply not implimented ",-2);
385 // return(-1); // unreachable
386}
387
388//=========================================================================
389
391
392 Epetra_Util util;
393 double* arrayPtr = A_;
394 const int size = (N_ == 1)? 1 : 4*(N_-1);
395 for(int j = 0; j < size ; j++) {
396 *arrayPtr++ = util.RandomDouble();
397 }
398
399 return(0);
400}
401
402void Ifpack_SerialTriDiMatrix::Print(std::ostream& os) const {
403 os <<" square format:"<<std::endl;
404 if(! A_ ) {
405 os <<" empty matrix "<<std::endl;
406 return;
407 }
408 for(int i=0 ; i < N_ ; ++i) {
409 for(int j=0 ; j < N_ ; ++j) {
410 if ( j >= i-1 && j <= i+1) {
411 os << (*this)(i,j)<<" ";
412 }
413 else
414 os <<". ";
415 }
416 os << std::endl;
417 }
418 }
Epetra_DataAccess
View
Copy
void UpdateFlops(int Flops_in) const
virtual void SetLabel(const char *const Label)
virtual int ReportError(const std::string Message, int ErrorCode) const
double RandomDouble()
Ifpack_SerialTriDiMatrix: A class for constructing and using real double precision general TriDi matr...
virtual const char * Label() const
Returns a character string describing the operator.
virtual double NormOne() const
Computes the 1-Norm of the this matrix.
int Shape(int NumRowCol)
Set dimensions of a Ifpack_SerialTriDiMatrix object; init values to zero.
int N() const
Returns column dimension of system.
double * A() const
Returns pointer to the this matrix.
virtual ~Ifpack_SerialTriDiMatrix()
Ifpack_SerialTriDiMatrix destructor.
int Random()
Column access function.
int Scale(double ScalarA)
Matrix-Vector multiplication, y = A*x, where 'this' == A.
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
bool operator==(const Ifpack_SerialTriDiMatrix &rhs) const
Comparison operator.
Ifpack_SerialTriDiMatrix & operator=(const Ifpack_SerialTriDiMatrix &Source)
Value copy from one matrix to another.
int Multiply(char TransA, char TransB, double ScalarAB, const Ifpack_SerialTriDiMatrix &A, const Ifpack_SerialTriDiMatrix &B, double ScalarThis)
Matrix-Matrix multiplication, this = ScalarThis*this + ScalarAB*A*B.
Ifpack_SerialTriDiMatrix(bool set_object_label=true)
Default constructor; defines a zero size object.
virtual double NormInf() const
Computes the Infinity-Norm of the this matrix.
Ifpack_SerialTriDiMatrix & operator+=(const Ifpack_SerialTriDiMatrix &Source)
Add one matrix to another.