00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "LinAlg.h"
00030 #include <math.h>
00031
00032 namespace linalg
00033 {
00034 using namespace linalg;
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 class MatrixPivoting : public Matrix
00045 {
00046 typedef REAL * INDEX;
00047 INDEX * const row_index;
00048
00049
00050 INDEX * const col_index;
00051
00052
00053
00054
00055
00056 double pivot_value;
00057 INDEX pivot_row;
00058 INDEX pivot_col;
00059 int pivot_odd;
00060
00061
00062 void pick_up_pivot(void);
00063
00064
00065 MatrixPivoting(const MatrixPivoting&);
00066 void operator = (const MatrixPivoting&);
00067
00068 public:
00069 MatrixPivoting(const Matrix& m);
00070
00071 ~MatrixPivoting(void);
00072
00073 double pivoting_and_elimination(void);
00074
00075
00076 };
00077
00078
00079
00080
00081
00082
00083
00084 MatrixPivoting::MatrixPivoting(const Matrix& m)
00085 : Matrix(m), row_index(new INDEX[nrows]),
00086 col_index(new INDEX[ncols]), pivot_value(0),
00087 pivot_row(0), pivot_col(0), pivot_odd(false)
00088 {
00089 assert( row_index != 0 && col_index != 0);
00090 register INDEX rp = elements;
00091 for(register INDEX * rip = row_index; rip<row_index+nrows;)
00092 *rip++ = rp++;
00093 register INDEX cp = elements;
00094 for(register INDEX * cip = col_index; cip<col_index+ncols; cp += nrows)
00095 *cip++ = cp;
00096 }
00097
00098 MatrixPivoting::~MatrixPivoting(void)
00099 {
00100 is_valid();
00101 delete row_index;
00102 delete col_index;
00103 }
00104
00105
00106
00107
00108
00109
00110
00111
00112 void MatrixPivoting::pick_up_pivot(void)
00113 {
00114 register REAL max_elem = -1;
00115 INDEX * prpp = 0;
00116 INDEX * pcpp = 0;
00117
00118 int col_odd = 0;
00119
00120 for(register INDEX * cpp = col_index; cpp < col_index + ncols; cpp++)
00121 {
00122 register const REAL * cp = *cpp;
00123 if( cp == 0 )
00124 continue;
00125 int row_odd = 0;
00126 for(register INDEX * rip = row_index; rip < row_index + nrows; rip++,cp++)
00127 if( *rip )
00128 {
00129 const REAL v = *cp;
00130 if( fabs(v) > max_elem )
00131 {
00132 max_elem = fabs(v);
00133 pivot_value = v;
00134 prpp = rip;
00135 pcpp = cpp;
00136 pivot_odd = row_odd ^ col_odd;
00137 }
00138 row_odd ^= 1;
00139 }
00140 col_odd ^= 1;
00141 }
00142
00143 assure( max_elem >= 0 && prpp != 0 && pcpp != 0,
00144 "All the rows and columns have been already pivoted and eliminated");
00145
00146
00147 pivot_row = *prpp, *prpp = 0;
00148 pivot_col = *pcpp, *pcpp = 0;
00149 }
00150
00151
00152
00153
00154
00155
00156 double MatrixPivoting::pivoting_and_elimination(void)
00157 {
00158 pick_up_pivot();
00159 if( pivot_value == 0 )
00160 return 0;
00161
00162 assert( pivot_row != 0 && pivot_col != 0 );
00163
00164 register REAL * pcp;
00165 register const INDEX * rip;
00166
00167 for(pcp=pivot_col,rip=row_index; rip<row_index+nrows; pcp++,rip++)
00168 if( *rip )
00169 *pcp /= pivot_value;
00170
00171
00172
00173
00174 const REAL * prp = pivot_row;
00175 for(register const INDEX * cpp = col_index; cpp < col_index + ncols; cpp++,prp+=nrows)
00176 {
00177 register REAL * cp = *cpp;
00178 if( cp == 0 )
00179 continue;
00180 const double fac = *prp;
00181
00182 for(pcp=pivot_col,rip=row_index; rip<row_index+nrows; pcp++,cp++,rip++)
00183 if( *rip )
00184 *cp -= *pcp * fac;
00185 }
00186
00187 return pivot_odd ? -pivot_value : pivot_value;
00188 }
00189
00190
00191
00192
00193
00194
00195
00196 double Matrix::determinant(void) const
00197 {
00198 is_valid();
00199
00200 if( nrows != ncols )
00201 info(), _error("Can't obtain determinant of a non-square matrix");
00202
00203 if( row_lwb != col_lwb )
00204 info(), _error("Row and col lower bounds are inconsistent");
00205
00206 MatrixPivoting mp(*this);
00207
00208 register double det = 1;
00209 register int k;
00210
00211 for(k=0; k<ncols && det != 0; k++)
00212 det *= mp.pivoting_and_elimination();
00213
00214 return det;
00215 }
00216 }