Ifpack Package Browser (Single Doxygen Collection) Development
ifp_DenseMat.h
Go to the documentation of this file.
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 #ifndef _IFP_DENSEMAT_H_
00031 #define _IFP_DENSEMAT_H_
00032 
00033 #include "Ifpack_config.h"
00034 
00035 #include <math.h>
00036 #include <assert.h>
00037 #include "ifp_LocalMat.h"
00038 #include "ifp_BlockVec.h"
00039 #include "ifp_LocalPrecon.h"
00040 #include "ifp_ifpack.h"
00041 #include "Epetra_BLAS.h"
00042 #include "Epetra_LAPACK.h"
00043 
00044 // ifp_DenseMat objects typically allocate their own storage for matrix data.
00045 // The default constructor is the only constructor that does not allocate
00046 // storage.  Thus this constructor may be used to control memory
00047 // allocation, for example by allocating an array of ifp_DenseMat using new.
00048 // In this case, the data pointer must be set to NULL before the destructor
00049 // is called.
00050 
00051 class IFPACK_DEPRECATED ifp_DenseMat : public ifp_LocalMat
00052 {
00053 // Temporary change protected:
00054 public:
00055     double *a;
00056     int    nrow;
00057     int    ncol;
00058 
00059 public:
00060     ifp_DenseMat() {a = NULL; nrow = ncol = 0;}
00061    ~ifp_DenseMat() {delete [] a;}
00062 
00063     ifp_DenseMat(const int r, const int c)
00064         {nrow = r; ncol = c; a = new double[r*c];}
00065 
00066     ifp_DenseMat(const ifp_DenseMat& A);
00067 // cannot inline on Cray
00068 //      {nrow = A.nrow; 
00069 //       ncol = A.ncol;
00070 //       register double *p = a = new double[nrow*ncol];
00071 //       register double *q = A.a;
00072 //       for (int i=0; i<nrow*ncol; i++) *p++ = *q++;}
00073 
00074     void set(const int r, const int c, double *d)
00075         {nrow = r; ncol = c; a = d;}
00076 
00077     int numrow() const {return nrow;}
00078     int numcol() const {return ncol;}
00079 
00080     double operator=(double s)
00081         {register double *p = a;
00082          for (int i=0; i<nrow*ncol; i++) *p++ = s; return s;}
00083 
00084     // 0-based indexing
00085     const double& operator()(const unsigned int i, const unsigned int j) const
00086         {return a[j*nrow+i];}
00087     double& operator()(const unsigned int i, const unsigned int j)
00088         {return a[j*nrow+i];}
00089 
00090     // virtual functions
00091 
00092     double *& Data() {return a;}
00093     const double *Data() const {return a;}
00094     ifp_LocalMat *CreateEmpty() const {return new ifp_DenseMat();}
00095     inline ifp_LocalMat *CreateInv(ifp_LocalPrecon&) const;
00096     inline void SetToZero(int, int);
00097     void MatCopy(const ifp_LocalMat& A)  // UNDONE: realloc if necessary
00098         {register double *p = a;
00099          register double *q = ((ifp_DenseMat *)&A)->a;
00100          for (int i=0; i<nrow*ncol; i++) *p++ = *q++;}
00101     void Print(std::ostream&) const;
00102 
00103     // virtual mathematical functions
00104 
00105     inline void Mat_Trans(ifp_LocalMat *B) const;
00106     inline void Mat_Mat_Add(const ifp_LocalMat *B, ifp_LocalMat *C, double alpha) const;
00107     inline void Mat_Mat_Mult(const ifp_LocalMat *B, ifp_LocalMat *C, 
00108         double alpha, double beta) const;
00109     inline void Mat_Vec_Mult(const ifp_BlockVec& B, ifp_BlockVec& C, 
00110         double alpha, double beta) const;
00111     inline void Mat_Trans_Vec_Mult(const ifp_BlockVec& B, ifp_BlockVec& C,
00112         double alpha, double beta) const;
00113     inline void Mat_Vec_Solve(const ifp_BlockVec& B, ifp_BlockVec& C) const;
00114     inline void Mat_Trans_Vec_Solve(const ifp_BlockVec& B, ifp_BlockVec& C) const;
00115 };
00116 
00117 std::ostream& operator << (std::ostream& os, const ifp_DenseMat& mat);
00118 
00119 class IFPACK_DEPRECATED ifp_DenseMat_LU : public ifp_DenseMat
00120 {
00121 private:
00122     double  *lu;
00123     integer *ipiv;  // FORTRAN integer
00124 
00125 public:
00126     inline ifp_DenseMat_LU(const ifp_DenseMat&);
00127     inline void Mat_Vec_Solve(const ifp_BlockVec&, ifp_BlockVec&) const;
00128     inline ~ifp_DenseMat_LU();
00129 };
00130 
00131 class IFPACK_DEPRECATED ifp_DenseMat_INVERSE : public ifp_DenseMat
00132 {
00133 private:
00134 public:
00135     inline ifp_DenseMat_INVERSE(const ifp_DenseMat&);
00136     void Mat_Vec_Solve(const ifp_BlockVec& B, ifp_BlockVec& X) const
00137         {solve_is_mult(B, X);}
00138 };
00139 
00140 class IFPACK_DEPRECATED ifp_DenseMat_SVD : public ifp_DenseMat
00141 {
00142 private:
00143 public:
00144     inline ifp_DenseMat_SVD(const ifp_DenseMat&, double, double);
00145     void Mat_Vec_Solve(const ifp_BlockVec& B, ifp_BlockVec& X) const
00146         {solve_is_mult(B, X);}
00147 };
00148 
00149 class IFPACK_DEPRECATED ifp_DenseMat_DIAGDOM : public ifp_DenseMat
00150 {
00151 private:
00152 public:
00153     inline ifp_DenseMat_DIAGDOM(const ifp_DenseMat&, double alpha);
00154     void Mat_Vec_Solve(const ifp_BlockVec& B, ifp_BlockVec& X) const
00155         {solve_is_mult(B, X);}
00156 };
00157 
00158 class IFPACK_DEPRECATED ifp_DenseMat_GERSH : public ifp_DenseMat
00159 {
00160 private:
00161 public:
00162     inline ifp_DenseMat_GERSH(const ifp_DenseMat&, double);
00163     void Mat_Vec_Solve(const ifp_BlockVec& B, ifp_BlockVec& X) const
00164         {solve_is_mult(B, X);}
00165 };
00166 
00167 // B = A'
00168 // must copy data
00169 inline IFPACK_DEPRECATED void ifp_DenseMat::Mat_Trans(ifp_LocalMat *B) const
00170 {
00171     ifp_DenseMat& b = *(ifp_DenseMat *) B;
00172 
00173     assert (this != &b);       // in-place not allowed
00174 
00175     if (b.a == NULL && b.nrow == 0 && b.ncol == 0)
00176     {
00177   b.nrow = ncol;
00178   b.ncol = nrow;
00179   b.a = new double[b.nrow*b.ncol];
00180     }
00181 
00182     assert (a != NULL);
00183     assert (b.a != NULL);
00184     assert (b.nrow == ncol);
00185     assert (b.ncol == nrow);
00186 
00187     int i, j;
00188     double *p, *q;
00189     
00190     // Traverse B and fill from A.
00191 
00192     p = b.a;
00193     for (i=0; i<nrow; i++)
00194     {
00195   q = a+i;
00196         for (j=0; j<ncol; j++)
00197   {
00198       *p++ = *q;
00199       q += nrow;
00200   }
00201     }
00202 }
00203 
00204 // C = A + alpha B
00205 inline IFPACK_DEPRECATED void ifp_DenseMat::Mat_Mat_Add(const ifp_LocalMat *B, ifp_LocalMat *C,
00206     double alpha) const
00207 {
00208     ifp_DenseMat& b = *(ifp_DenseMat *) B;
00209     ifp_DenseMat& c = *(ifp_DenseMat *) C;
00210 
00211     if (c.a == NULL && c.nrow == 0 && c.ncol == 0)
00212     {
00213   c.nrow = nrow;
00214   c.ncol = ncol;
00215   c.a = new double[c.nrow*c.ncol];
00216     }
00217 
00218     assert (a != NULL);
00219     assert (b.a != NULL);
00220     assert (c.a != NULL);
00221     assert (nrow == b.nrow);
00222     assert (ncol == b.ncol);
00223     assert (nrow == c.nrow);
00224     assert (ncol == c.ncol);
00225 
00226     int i;
00227     double *ap = a;
00228     double *bp = b.a;
00229     double *cp = c.a;
00230 
00231     if (alpha == 1.0)
00232         for (i=0; i<nrow*ncol; i++)
00233             *cp++ = *ap++ + *bp++;
00234 
00235     else if (alpha == -1.0)
00236         for (i=0; i<nrow*ncol; i++)
00237       *cp++ = *ap++ - *bp++;
00238 
00239     else
00240         for (i=0; i<nrow*ncol; i++)
00241       *cp++ = *ap++ + alpha * *bp++;
00242 }
00243 
00244 // C = alpha A B + beta C
00245 inline IFPACK_DEPRECATED void ifp_DenseMat::Mat_Mat_Mult(const ifp_LocalMat *B, ifp_LocalMat *C, 
00246     double alpha, double beta) const
00247 {
00248     assert (B != C);
00249     ifp_DenseMat& b = *(ifp_DenseMat *) B;
00250     ifp_DenseMat& c = *(ifp_DenseMat *) C;
00251 
00252     if (beta == 0.0 && c.a == NULL && c.nrow == 0 && c.ncol == 0)
00253     {
00254   c.nrow = nrow;
00255   c.ncol = b.ncol;
00256   c.a = new double[c.nrow*c.ncol];
00257     }
00258 
00259     if (beta == 0.0 && (c.nrow != nrow || c.ncol != b.ncol))
00260     {
00261         assert (c.a != NULL);
00262   delete c.a;
00263 
00264   c.nrow = nrow;
00265   c.ncol = b.ncol;
00266   c.a = new double[c.nrow*c.ncol];
00267     }
00268 
00269     assert (a != NULL);
00270     assert (b.a != NULL);
00271     assert (c.a != NULL);
00272 
00273     char transa = 'N';
00274     char transb = 'N';
00275     integer M = nrow;
00276     integer N = c.ncol;
00277     integer K = ncol;
00278     integer LDA = nrow;
00279     integer LDB = b.nrow;
00280     integer LDC = c.nrow;
00281 
00282     Epetra_BLAS blas;
00283     blas.GEMM(transa, transb, M, N, K, alpha, a, LDA, b.a, LDB, beta, c.a, LDC);
00284     //F77NAME(dgemm)(&transa, &transb, &M, &N, &K, &alpha, a, &LDA, 
00285     //    b.a, &LDB, &beta, c.a, &LDC);
00286 }
00287 
00288 // C = alpha A B + beta C
00289 inline IFPACK_DEPRECATED void ifp_DenseMat::Mat_Vec_Mult(const ifp_BlockVec& B, ifp_BlockVec& C,
00290     double alpha, double beta) const
00291 {
00292     char trans = 'N';
00293     integer M = nrow;
00294     integer N = C.dim1;
00295     integer K = ncol;
00296     integer LDA = nrow;
00297     integer LDB = B.dim0;
00298     integer LDC = C.dim0;
00299 
00300     assert (a != NULL);
00301     assert (&B != &C);
00302 
00303     Epetra_BLAS blas;
00304     blas.GEMM(trans, trans, M, N, K, alpha, a, LDA, B.v, LDB, beta, C.v, LDC);
00305     //F77NAME(dgemm)(&trans, &trans, &M, &N, &K, &alpha, a, &LDA, 
00306     //    B.v, &LDB, &beta, C.v, &LDC);
00307 }
00308 
00309 inline IFPACK_DEPRECATED void ifp_DenseMat::Mat_Trans_Vec_Mult(const ifp_BlockVec& B, 
00310     ifp_BlockVec& C, double alpha, double beta) const
00311 {
00312     char transa = 'T';
00313     char transb = 'N';
00314     integer M = ncol;
00315     integer N = C.dim1;
00316     integer K = nrow;
00317     integer LDA = nrow;
00318     integer LDB = B.dim0;
00319     integer LDC = C.dim0;
00320 
00321     assert (a != NULL);
00322     assert (&B != &C);
00323 
00324     Epetra_BLAS blas;
00325     blas.GEMM(transa, transb, M, N, K, alpha, a, LDA, B.v, LDB, beta, C.v, LDC);
00326     //F77NAME(dgemm)(&transa, &transb, &M, &N, &K, &alpha, a, &LDA,
00327     //    B.v, &LDB, &beta, C.v, &LDC);
00328 }
00329 
00330 inline IFPACK_DEPRECATED void ifp_DenseMat::Mat_Vec_Solve(const ifp_BlockVec&, ifp_BlockVec&) const
00331 {
00332     ifp_error("ifp_DenseMat::Mat_Vec_Solve: to solve with a local matrix,\n"
00333             "a local preconditioner should be used.\n", 0);
00334 }
00335 
00336 inline IFPACK_DEPRECATED void ifp_DenseMat::Mat_Trans_Vec_Solve(const ifp_BlockVec&,
00337 ifp_BlockVec&) const
00338 {
00339     ifp_error("ifp_DenseMat::Mat_Trans_Vec_Solve: to solve with a local matrix,\n"
00340             "a local preconditioner should be used.\n", 0);
00341 }
00342 
00343 inline IFPACK_DEPRECATED void ifp_DenseMat::SetToZero(int r, int c)
00344 {
00345     if (a == NULL && nrow == 0 && ncol == 0)
00346     {
00347   nrow = r;
00348   ncol = c;
00349   a = new double[nrow*ncol];
00350     }
00351 
00352     assert (a != NULL);
00353     assert (nrow == r);
00354     assert (ncol == c);
00355 
00356     double *p = a;
00357     for (int i=0; i<nrow*ncol; i++)
00358   *p++ = 0.0;
00359 }
00360 
00361 inline IFPACK_DEPRECATED ifp_LocalMat *ifp_DenseMat::CreateInv(ifp_LocalPrecon& local_precon) const
00362 {
00363     assert (a != NULL);
00364 
00365     switch (local_precon.name)
00366     {
00367     case LP_LU:
00368         return new ifp_DenseMat_LU(*this);
00369     case LP_INVERSE:
00370         return new ifp_DenseMat_INVERSE(*this);
00371     case LP_SVD:
00372         return new ifp_DenseMat_SVD(*this, local_precon.darg1, local_precon.darg2);
00373     case LP_DIAGDOM:
00374         return new ifp_DenseMat_DIAGDOM(*this, local_precon.darg1);
00375     case LP_GERSH:
00376         return new ifp_DenseMat_GERSH(*this, local_precon.darg1);
00377     default:
00378         ifp_error("The local preconditioner you have chosen is not available\n"
00379                 "for the type of blocks (ifp_DenseMat) in the block matrix.", 0);
00380     }
00381     return(0); // Zero pointer (never returned because ifp_error aborts, but this satisfies the compiler)
00382 }
00383 
00384 inline IFPACK_DEPRECATED ifp_DenseMat_LU::ifp_DenseMat_LU(const ifp_DenseMat& A)
00385 {
00386     a = NULL; // indicate this is implied inverse
00387     nrow = A.numrow();
00388     ncol = A.numcol();
00389 
00390     assert (nrow == ncol);
00391 
00392     lu = new double[nrow*ncol];
00393     // copy to lu
00394 
00395     double *p = lu;
00396     const double *q = &A(0,0);
00397     for (int i=0; i<nrow*ncol; i++)
00398         *p++ = *q++;
00399 
00400     integer M = nrow;
00401     integer N = ncol;
00402     integer LDA = nrow;
00403     integer INFO;
00404 
00405     ipiv = new integer[nrow];
00406 
00407     Epetra_LAPACK lapack;
00408     lapack.GETRF(M, N, lu, LDA, ipiv, &INFO);
00409     //F77NAME(dgetrf)(&M, &N, lu, &LDA, ipiv, &INFO);
00410     if (INFO != 0)
00411   ifp_error("ifp_DenseMat_LU: dgetrf error", INFO);
00412 }
00413 
00414 inline IFPACK_DEPRECATED void ifp_DenseMat_LU::Mat_Vec_Solve(const ifp_BlockVec& b, ifp_BlockVec& x) const
00415 {
00416     char trans = 'N';
00417     integer N = ncol;
00418     integer NRHS = b.dim1;
00419     integer LDA = nrow;
00420     integer LDB = b.dim0;
00421     integer INFO;
00422 
00423     if (b.v != x.v)
00424         x.BlockCopy(b);
00425 
00426     Epetra_LAPACK lapack;
00427     lapack.GETRS(trans, N, NRHS, lu, LDA, ipiv, x.v, LDB, &INFO);
00428     //F77NAME(dgetrs)(&trans, &N, &NRHS, lu, &LDA, ipiv, 
00429     //    x.v, &LDB, &INFO);
00430     if (INFO != 0)
00431   ifp_error("ifp_DenseMat_LU: dgetrs error", INFO);
00432 }
00433 
00434 inline IFPACK_DEPRECATED ifp_DenseMat_LU::~ifp_DenseMat_LU()
00435 {
00436     delete [] ipiv;
00437     delete [] lu;
00438 }
00439 
00440 inline IFPACK_DEPRECATED ifp_DenseMat_INVERSE::ifp_DenseMat_INVERSE(const ifp_DenseMat& A)
00441     : ifp_DenseMat(A)
00442 {
00443     integer M = nrow;
00444     integer N = ncol;
00445     integer LDA = nrow;
00446     integer INFO;
00447 
00448     assert (nrow == ncol);
00449 
00450     integer *ipiv = new integer[nrow];
00451     integer LWORK = N*N;
00452     double *work = new double[LWORK];
00453 
00454     // LU factorize
00455     Epetra_LAPACK lapack;
00456     lapack.GETRF(M, N, a, LDA, ipiv, &INFO);
00457     //F77NAME(dgetrf)(&M, &N, a, &LDA, ipiv, &INFO);
00458     if (INFO != 0)
00459   ifp_error("ifp_DenseMat_INVERSE: dgetrf error", INFO);
00460 
00461     // compute inverse
00462     lapack.GETRI(N, a, LDA, ipiv, work, &LWORK, &INFO);
00463     //F77NAME(dgetri)(&N, a, &LDA, ipiv, work, &LWORK, &INFO);
00464     if (INFO != 0)
00465   ifp_error("ifp_DenseMat_INVERSE: dgetri error", INFO);
00466 
00467     delete [] ipiv;
00468     delete [] work;
00469 }
00470 
00471 inline IFPACK_DEPRECATED ifp_DenseMat_SVD::ifp_DenseMat_SVD(const ifp_DenseMat& A, double rthresh, 
00472   double athresh)
00473     : ifp_DenseMat(A)
00474 {
00475     assert (nrow == ncol);
00476 
00477     int i, j;
00478     double thresh;
00479     double *u  = new double[nrow*ncol];
00480     double *s  = new double[nrow];
00481     double *vt = new double[nrow*ncol];
00482 
00483     char job = 'A';
00484     integer N = nrow;
00485     integer LWORK = 5*N;
00486     double *work = new double[LWORK];
00487     integer INFO;
00488     integer nzero_diag = 0;
00489     for (i=0;i<nrow; i++) if (a[i*N+i] == 0.0) nzero_diag++;
00490 
00491      Epetra_LAPACK lapack;
00492     lapack.GESVD(job, job, N, N, a, N, s, u, N, vt, N, work, &LWORK, &INFO);
00493     //F77NAME(dgesvd) (&job, &job, &N, &N, a, &N, s, u, &N, vt, &N,
00494     //        work, &LWORK, &INFO);
00495     //if (INFO != 0)
00496     if (INFO < 0)
00497   ifp_error("ifp_DenseMat_SVD: dgesvd error", INFO);
00498 
00499     delete [] work;
00500 
00501     // apply threshold
00502     thresh = s[0]*rthresh + athresh;
00503     int n_replaced = 0;
00504     for (i=0; i<nrow; i++)
00505   if (s[i] < thresh) 
00506     {
00507       s[i] = thresh;
00508       n_replaced++;
00509     }
00510     //   cout <<nzero_diag<<" "<<s0<<" "<<sn<<" "<<s0/sn<<" "<<thresh<<" "<<n_replaced << endl;
00511 
00512 
00513     if (s[0] == 0.0)
00514   ifp_error("ifp_DenseMat_SVD: block is zero after thresholding", 0);
00515 
00516     // scale the columns of u with reciprocal singular values
00517     double *p = u;
00518     for (i=0; i<ncol; i++)
00519     {
00520   double scal = s[i];
00521   for (j=0; j<nrow; j++)
00522       *p++ /= scal;
00523     }
00524 
00525     // multiply through and store the result
00526     char trans = 'T';
00527     double alpha = 1.0;
00528     double beta = 0.0;
00529     Epetra_BLAS blas;
00530     blas.GEMM(trans, trans, N, N, N, alpha, vt, N, u, N, beta, a, N);
00531     //F77NAME(dgemm)(&trans, &trans, &N, &N, &N, &alpha, vt, &N,
00532     //    u, &N, &beta, a, &N);
00533 
00534     delete [] u;
00535     delete [] s;
00536     delete [] vt;
00537 }
00538 
00539 inline IFPACK_DEPRECATED ifp_DenseMat_DIAGDOM::ifp_DenseMat_DIAGDOM(const ifp_DenseMat& A, double alpha)
00540     : ifp_DenseMat(A)
00541 {
00542     int i, j;
00543     integer M = nrow;
00544     integer N = ncol;
00545     integer LDA = nrow;
00546     integer INFO;
00547 
00548     assert (nrow == ncol);
00549 
00550     integer *ipiv = new integer[nrow];
00551     integer LWORK = N*N;
00552     double *work = new double[LWORK];
00553 
00554     // compute sum of abs of off-diagonals and put in work array
00555     double *p = a;
00556     for (i=0; i<nrow; i++)
00557         work[i] = 0.0;
00558 
00559     for (j=0; j<ncol; j++)
00560     {
00561         for (i=0; i<nrow; i++)
00562         {
00563             if (i != j)
00564             {
00565                 work[i] += ABS(*p);
00566             }
00567             p++;
00568         }
00569     }
00570 
00571     for (i=0; i<nrow; i++)
00572     {
00573   if (ABS(a[i*nrow+i]) < alpha*work[i])
00574       a[i*nrow+i] = SGN(a[i*nrow+i])*alpha*work[i];
00575     }
00576 
00577     // LU factorize
00578     Epetra_LAPACK lapack;
00579     lapack.GETRF(M, N, a, LDA, ipiv, &INFO);
00580     //F77NAME(dgetrf)(&M, &N, a, &LDA, ipiv, &INFO);
00581     if (INFO != 0)
00582   ifp_error("ifp_DenseMat_DIAGDOM: dgetrf error", INFO);
00583 
00584     // compute inverse
00585     lapack.GETRI(N, a, LDA, ipiv, work, &LWORK, &INFO);
00586     //F77NAME(dgetri)(&N, a, &LDA, ipiv, work, &LWORK, &INFO);
00587     if (INFO != 0)
00588   ifp_error("ifp_DenseMat_DIAGDOM: dgetri error", INFO);
00589 
00590     delete [] ipiv;
00591     delete [] work;
00592 }
00593 
00594 inline IFPACK_DEPRECATED ifp_DenseMat_GERSH::ifp_DenseMat_GERSH(const ifp_DenseMat& A, double alpha)
00595     : ifp_DenseMat(A)
00596 {
00597     int i, j;
00598     integer M = nrow;
00599     integer N = ncol;
00600     integer LDA = nrow;
00601     integer INFO;
00602 
00603     assert (nrow == ncol);
00604 
00605     integer *ipiv = new integer[nrow];
00606     integer LWORK = N*N;
00607     double *work = new double[LWORK];
00608 
00609     // compute sum of abs of off-diagonals and put in work array
00610     double *p = a;
00611     for (i=0; i<nrow; i++)
00612         work[i] = 0.0;
00613 
00614     for (j=0; j<ncol; j++)
00615     {
00616         for (i=0; i<nrow; i++)
00617         {
00618             if (i != j)
00619             {
00620                 work[i] += ABS(*p);
00621             }
00622             p++;
00623         }
00624     }
00625 
00626     double aii;
00627     for (i=0; i<nrow; i++)
00628     {
00629   aii = a[i*nrow+i];
00630 
00631   if (aii >= 0.0)
00632         {
00633       if (aii - work[i] < alpha)
00634     a[i*nrow+i] += alpha - aii + work[i];
00635   }
00636   else
00637   {
00638       if (aii + work[i] > -alpha)
00639     a[i*nrow+i] -= alpha + aii + work[i];
00640   }
00641     }
00642 
00643     // LU factorize
00644     Epetra_LAPACK lapack;
00645     lapack.GETRF(M, N, a, LDA, ipiv, &INFO);
00646     //F77NAME(dgetrf)(&M, &N, a, &LDA, ipiv, &INFO);
00647     if (INFO != 0)
00648   ifp_error("ifp_DenseMat_GERSH: dgetrf error", INFO);
00649 
00650     // compute inverse
00651     lapack.GETRI(N, a, LDA, ipiv, work, &LWORK, &INFO);
00652     //F77NAME(dgetri)(&N, a, &LDA, ipiv, work, &LWORK, &INFO);
00653     if (INFO != 0)
00654   ifp_error("ifp_DenseMat_GERSH: dgetri error", INFO);
00655 
00656     delete [] ipiv;
00657     delete [] work;
00658 }
00659 
00660 #endif // _IFP_DENSEMAT_H_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines