00001 /***************************************************************************
00002                           lusolve.h  -  description
00003                              -------------------
00004     begin                : Mon May 8 2000
00005     copyright            : (C) 2000 by Michael Peeters
00006     email                : Michael.Peeters@vub.ac.be
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #ifndef LUSOLVE_H
00019 #define LUSOLVE_H
00020 
00021 #include "utility.h"
00022 #include <stdexcept>
00023 #include "numerictraits.h"
00024 #include "numerictypes.h"
00025 
00026 namespace MODEL {
00027 
00028   template <integer dims, class NT = NumericTraits<number,dims> >
00029   class LUSolve
00030   {
00031   public:
00032         typedef typename NT::number     numT;
00033         typedef typename NT::vect       vect;   
00034         typedef typename NT::matrix     matrix;
00035                 
00036         LUSolve(matrix& m, bool own=false) : M(own?m.clone():&m), owned(own) {decompose();};
00037         ~LUSolve() {if(owned) delete M;}
00038                 
00040         vect    operator()(const vect& u) {vect res(u); solve(res); return
00041                                                                                                                                   res;}
00042 
00044         void    solve(vect& u);
00046         numT    det(void);
00047                 
00048   private:
00049         // to avoid silly errors
00050         NO_COPY(LUSolve);
00051                 
00052         // perform decomposition
00053         void    decompose(void);
00054 
00055   private:
00056         matrix*                                 M;              
00057         integer                                 pivotrows[dims];
00058         numT                                    pivotsign;
00059                         
00060         bool    owned;
00061 
00062         static const numT       tiny=1.E-20;
00063   };
00064 
00065   template <integer dims, class NT >
00066   void LUSolve<dims,NT>::decompose(void)
00067   {
00068     matrix&             a(*M);
00069 
00070         pivotsign=1.0;  
00071         vect    rowscale(1.);
00072         // Test for singularity
00073         for (integer i=0;i<dims;i++)
00074           {
00075                 numT big=0.0;
00076                 numT t(0.);
00077                 for (integer j=0;j<dims;j++)
00078                   if ((t=fabs(a[i][j])) > big) big=t;
00079                 if (big == 0.0) throw std::logic_error("Singular matrix in routine ludcmp");
00080                 rowscale[i]/=big;
00081           }
00082         
00083         numT sum(0.);
00084         for (integer j=0;j<dims;j++)
00085           {
00086                 for (integer i=0;i<j;i++)
00087                   {
00088                         sum=a[i][j];
00089                         for (integer k=0;k<i;k++) sum -= a[i][k]*a[k][j];
00090                         a[i][j]=sum;
00091                   }
00092                 numT    big(0.);
00093                 integer imax(-1);
00094                 numT    dum(0.);
00095                 for (integer i=j;i<dims;i++)
00096                   {
00097                         sum=a[i][j];
00098                         for (integer k=0;k<j;k++) sum -= a[i][k]*a[k][j];
00099                         a[i][j]=sum;
00100                         if ( (dum=rowscale[i]*abs(sum) ) >= big) {
00101                           big=dum;
00102                           imax=i;
00103                         }
00104                   }
00105                 if (j != imax) {
00106                   for (integer k=0;k<dims;k++) {
00107                         dum=a[imax][k];
00108                         a[imax][k]=a[j][k];
00109                         a[j][k]=dum;
00110                   }
00111                   pivotsign *= -1;
00112                   rowscale[imax]=rowscale[j];
00113                 }
00114                 pivotrows[j]=imax;
00115                 if (a[j][j] == 0.0) a[j][j]=tiny;
00116                 if (j != (dims-1)) {
00117                   dum=1.0/(a[j][j]);
00118                   for (integer i=j+1;i<dims;i++) a[i][j] *= dum;
00119                 }
00120           }
00121   }
00122 
00123   template <integer dims, class NT>
00124   void LUSolve<dims,NT>::solve(vect& u)
00125   {
00126     matrix&             a(*M);
00127     numT                sum(0.);
00128         integer         ii(-1);
00129         bool    nonzero=false;
00130         
00131         for (integer i=0;i<dims;i++)
00132           {             
00133                 integer ip=pivotrows[i];
00134                 sum=u[ip];
00135                 u[ip]=u[i];
00136                 
00137                 if (nonzero) for (integer j=ii;j<i;j++) sum -= a[i][j]*u[j];
00138                 else if (sum!=0.0) {nonzero=true; ii=i;}
00139                 
00140                 u[i]=sum;
00141           }
00142         for (integer i(dims-1);i>=0;i--) {
00143           sum = u[i];
00144           for (integer j(i+1);j<dims;j++) sum -= a[i][j]*u[j];
00145           u[i]=sum/a[i][i];
00146         }
00147   }
00148 
00149   template <integer dims, class NT>
00150   LUSolve<dims,NT>::numT        
00151   LUSolve<dims,NT>::det(void)
00152   {
00153         numT    d(pivotsign);
00154         for(integer j=0;j<dims;j++) d *= (*M)[j][j];
00155         return d;
00156   }
00157 
00158 
00159 } // end MODEL;
00160 
00161 #endif
00162 
00163 /*********************************************************************
00164 $Id: lusolve.h,v 1.1.2.1 2001/11/15 13:33:26 mpeeters Exp $
00165 **********************************************************************
00166 
00167 $Log: lusolve.h,v $
00168 Revision 1.1.2.1  2001/11/15 13:33:26  mpeeters
00169 Last commit before moving everything to new gcc stdc++ headers (3.0).
00170 
00171 Revision 1.1  2001/05/22 10:54:55  mpeeters
00172 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter".
00173 
00174 Revision 1.4  2001/05/21 11:53:16  mpeeters
00175 Removed Makefile.in, which is automatically generated anyway.
00176 
00177 Revision 1.3  2000/09/22 08:54:08  mpeeters
00178 Changed dependencies. Missed this one on the last commit
00179 
00180 Revision 1.2  2000/09/15 10:26:31  mpeeters
00181 Cleaned out header and added CVS tails to files, removed superfuous
00182 @author comments, inserted dates
00183 
00184 *********************************************************************/