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