IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
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 //=============================================================================
147 int 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 //=============================================================================
209 void 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 //=============================================================================
313 void 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 //=============================================================================
364 int Ifpack_SerialTriDiMatrix::Scale(double ScalarA) {
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 // //=========================================================================
380 int 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 
402 void 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  }
int Shape(int NumRowCol)
Set dimensions of a Ifpack_SerialTriDiMatrix object; init values to zero.
bool operator==(const Ifpack_SerialTriDiMatrix &rhs) const
Comparison operator.
virtual double NormOne() const
Computes the 1-Norm of the this matrix.
double * A() const
Returns pointer to the this matrix.
Ifpack_SerialTriDiMatrix & operator+=(const Ifpack_SerialTriDiMatrix &Source)
Add one matrix to another.
virtual double NormInf() const
Computes the Infinity-Norm of the this matrix.
int N() const
Returns column dimension of system.
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream &lt;&lt; operator.
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.
double RandomDouble()
Ifpack_SerialTriDiMatrix & operator=(const Ifpack_SerialTriDiMatrix &Source)
Value copy from one matrix to another.
int Random()
Column access function.
int Scale(double ScalarA)
Matrix-Vector multiplication, y = A*x, where &#39;this&#39; == A.
Ifpack_SerialTriDiMatrix(bool set_object_label=true)
Default constructor; defines a zero size object.
Epetra_DataAccess
virtual ~Ifpack_SerialTriDiMatrix()
Ifpack_SerialTriDiMatrix destructor.
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.