ifp_BlockMat.cpp

00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "ifp_BlockMat.h"
00031 #include "ifp_BlockVec.h"
00032 #include "ifp_DenseMat.h"
00033 #include "ifp_ifpack.h"
00034 /* Epetra typedef no longer used */
00035 /*#include "Epetra_Object.h" // Bring in Epetra_fmtflags typedef*/
00036 /*#define DEBUG */
00037 
00038 ifp_BlockMat::~ifp_BlockMat()
00039 {
00040     if (a != NULL)
00041     {
00042       //delete [] a[0]->Data();       // values
00043             for (int i=0; i<nnz; i++)
00044                 a[i]->Data() = NULL;      // so DenseMat destructor safe
00045             delete [] (ifp_DenseMat *) a[0];  // must be array of ifp_DenseMat
00046     }
00047     delete [] a;
00048     //delete [] ja;
00049     //delete [] ia;
00050     //delete kvstr;
00051     //delete kvstc;
00052 }
00053 
00054 void ifp_BlockMat::mult(int nr, int nc, const double *B, int ldu, 
00055     double *C, int ldv) const
00056 {
00057     ifp_BlockVec Bb(nr, nc, B, ldu, kvstc);
00058     ifp_BlockVec Cb(dimrow(), nc, C, ldv, kvstr);
00059     Cb.VecSetToZero();
00060 
00061     for (int i=0; i<nrow; i++)
00062     {
00063         for (int j=ia[i]; j<ia[i+1]; j++)
00064     {
00065              a[j]->Mat_Vec_Mult(Bb(ja[j]), Cb(i), 1.0, 1.0);
00066     }
00067     }
00068 }
00069 
00070 void ifp_BlockMat::trans_mult(int nr, int nc, const double *B, int ldu,
00071     double *C, int ldv) const
00072 {
00073     ifp_BlockVec Bb(nr, nc, B, ldu, kvstc);
00074     ifp_BlockVec Cb(dimcol(), nc, C, ldv, kvstr);
00075     Cb.VecSetToZero();
00076 
00077     for (int i=0; i<nrow; i++)
00078     {
00079         for (int j=ia[i]; j<ia[i+1]; j++)
00080     {
00081              a[j]->Mat_Trans_Vec_Mult(Bb(i), Cb(ja[j]), 1.0, 1.0);
00082     }
00083     }
00084 }
00085 
00086 // Converts Aztec Matrix to IFPACK matrix
00087 
00088 ifp_BlockMat::ifp_BlockMat(double *val, int *indx, int * bindx, int *rpntr, int *cpntr,
00089            int *bpntr, int nrow_, int ncol_, int nnz_, int nnzs_)
00090 {
00091     int i, j, k;
00092     double *pp;
00093 
00094     a = (ifp_LocalMat **) NULL;
00095     ja = bindx;
00096     ia = bpntr;
00097     kvstr = rpntr;
00098     kvstc = cpntr;
00099     nrow = nrow_;
00100     ncol = ncol_;
00101     nnz = nnz_;
00102     nnzs = nnzs_;
00103 
00104 #ifdef DEBUG
00105      cerr << nrow << " block rows.\n";
00106      cerr << ncol << " block columns.\n";
00107 
00108      cerr << nnz << " block nonzeros.\n";
00109      cerr << nnzs << " stored scalar entries.\n";
00110 #endif
00111 
00112     a  = new ifp_LocalMatp[nnz];   // array of pointers to blocks
00113 
00114     ifp_DenseMat *p = new ifp_DenseMat[nnz]; // array of blocks
00115 
00116     for (k=0; k<nnz; k++)
00117         a[k] = &p[k];
00118 
00119     pp = val;  // array of values
00120 
00121     // fill entries and ja -----------------------------------------------------
00122 
00123     k = 0;
00124 
00125     // loop on block rows
00126     for (i=0; i<nrow; i++)
00127     {
00128         int neqr = kvstr[i+1]-kvstr[i];
00129 
00130     for (j=ia[i]; j<ia[i+1]; j++)
00131         {
00132       int neqc = kvstc[bindx[j]+1] - kvstc[bindx[j]];
00133       p[k].set(neqr,neqc,pp);
00134       /*      printf("In az2bp i = %d  j = %d\n",i,j);
00135           cout << p[k]; */
00136       pp = pp + neqr * neqc;
00137       k++;
00138         }
00139     }
00140 
00141 }
00142 
00143 // Non-member functions
00144 
00145 ostream& operator << (ostream& os, const ifp_BlockMat& mat)
00146 {
00147         int M = mat.numrow();
00148         int N = mat.numcol();
00149         int rowp1, colp1;
00150         int flag = 0;
00151 /*        Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
00152         Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
00153         int oldp = os.precision(12); */
00154 
00155         for (int i = 0; i < M ; i++)
00156            for (int j=mat.row_ptr(i);j<mat.row_ptr(i+1);j++)
00157            {   
00158               rowp1 =  i + 1;
00159               colp1 =  mat.col_ind(j) + 1;
00160               if ( rowp1 == M && colp1 == N ) flag = 1;
00161               os.width(14);
00162               os <<  rowp1 ; os << "    " ;
00163               os.width(14);
00164               os <<  colp1;
00165 #if 0
00166               os.width(20);
00167               os <<  *(mat.val(j).Data());
00168 #endif
00169               os << endl;
00170            }
00171 
00172 #if 0
00173         if (flag == 0)
00174         {
00175            os.width(14);
00176            os <<  M ; os << "    " ;
00177            os.width(14);
00178            os <<  N ; os << "    " ;
00179            os.width(20);
00180            os <<  mat(M-1,N-1) << "\n";
00181         }
00182 #endif
00183 
00184 /*        os.setf(olda,ios::adjustfield);
00185         os.setf(oldf,ios::floatfield);
00186         os.precision(oldp); */
00187 
00188         return os;
00189 }

Generated on Wed May 12 21:46:02 2010 for IFPACK by  doxygen 1.4.7