[6ead9d] | 1 | /* =================================================================== |
---|
| 2 | Dateiname: lgs.cc |
---|
| 3 | =================================================================== */ |
---|
| 4 | |
---|
| 5 | #include "lgs.h" |
---|
| 6 | #include "bifacConfig.h" |
---|
| 7 | |
---|
[0c4a34b] | 8 | #ifdef HAVE_BIFAC |
---|
[6ead9d] | 9 | |
---|
| 10 | //--<>--------------------------------- |
---|
| 11 | LGS::LGS( int r, int c, bool inv )// KONSTRUKTOR |
---|
| 12 | //--<>--------------------------------- |
---|
| 13 | { |
---|
| 14 | INVERSE = inv; |
---|
| 15 | |
---|
| 16 | max_columns = c; |
---|
[806c18] | 17 | max_rows = (r>c) ? c : r ; |
---|
[6ead9d] | 18 | |
---|
| 19 | if( INVERSE ) |
---|
| 20 | { |
---|
| 21 | c *=2; |
---|
| 22 | } |
---|
| 23 | |
---|
| 24 | CFMatrix AA(max_rows, c); |
---|
| 25 | CFMatrix bb(max_rows, 1); |
---|
| 26 | A=AA; |
---|
| 27 | b=bb; |
---|
| 28 | pivot = new int [max_columns+1]; // We start with index 1 (easier) |
---|
| 29 | for( int i=0; i<=max_columns; i++) pivot[i]=0; |
---|
| 30 | now_row = 1; // We start counting in the first row |
---|
| 31 | } |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | //--<>--------------------------------- |
---|
| 35 | LGS::~LGS( void )// DESTRUKTOR |
---|
| 36 | //--<>--------------------------------- |
---|
| 37 | { |
---|
| 38 | delete[] pivot; |
---|
[806c18] | 39 | } |
---|
[6ead9d] | 40 | |
---|
| 41 | //--<>--------------------------------- |
---|
| 42 | void LGS::reset(void) |
---|
| 43 | //--<>--------------------------------- |
---|
| 44 | { // Clear the matrix for a new computation |
---|
| 45 | now_row=1; |
---|
[806c18] | 46 | for( int i=0; i<=max_columns; i++) |
---|
[6ead9d] | 47 | pivot[i]=0; |
---|
| 48 | } |
---|
| 49 | |
---|
| 50 | //--<>--------------------------------- |
---|
| 51 | bool LGS::new_row( const CFMatrix Z, const CanonicalForm bb) |
---|
| 52 | //--<>--------------------------------- |
---|
[806c18] | 53 | { // Insert a new row |
---|
[6ead9d] | 54 | ASSERT ( (1 <= now_row && now_row <=max_rows), "wrong number of rows => Matrix has max. rank"); |
---|
| 55 | int i; |
---|
| 56 | |
---|
| 57 | // if (INVERSE) |
---|
[806c18] | 58 | // cout << "* Integriere Zeile "<<now_row << " (max " <<max_rows<<" x " |
---|
[6ead9d] | 59 | // << max_columns << ")\n" |
---|
| 60 | // << "Z = " << Z << "\nb = " << bb << endl << flush; |
---|
| 61 | |
---|
| 62 | |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | // === Insert a new row === |
---|
| 66 | for(i=1; i<=max_columns; i++) |
---|
| 67 | A(now_row, i) = Z(1,i); |
---|
[806c18] | 68 | b(now_row, 1) = bb; |
---|
| 69 | |
---|
[6ead9d] | 70 | |
---|
| 71 | // cout << "* reduzierte Matrix (vor lin_dep) " << A << endl; |
---|
| 72 | // === check linear dependency === |
---|
[806c18] | 73 | if ( ! lin_dep() ) |
---|
[6ead9d] | 74 | now_row++; |
---|
[806c18] | 75 | else |
---|
[6ead9d] | 76 | return(false); |
---|
[806c18] | 77 | |
---|
[6ead9d] | 78 | // === Reduce the previous rows === |
---|
| 79 | for(i=1; i<now_row-1; i++) |
---|
| 80 | reduce(now_row-1,i ); |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | // cout << "\n* Zeile Z =" << Z << " ist integriert!\n" << A << flush; |
---|
| 84 | |
---|
| 85 | // if( INVERSE ) cout << A; |
---|
| 86 | // cout << "* Verlasse new_row!\n" << "* Z = " << Z << endl << flush; |
---|
| 87 | return(true); // row was linear independent |
---|
[806c18] | 88 | } |
---|
| 89 | |
---|
[6ead9d] | 90 | |
---|
| 91 | //--<>--------------------------------- |
---|
| 92 | bool LGS::lin_dep( void ) |
---|
| 93 | //--<>--------------------------------- |
---|
| 94 | { // check if new row is linear dependent of the former ones |
---|
| 95 | |
---|
| 96 | int i; |
---|
| 97 | |
---|
| 98 | // === Reduce the actual row === |
---|
| 99 | for(i=1; i<now_row; i++) |
---|
| 100 | reduce(i, now_row); |
---|
| 101 | |
---|
| 102 | // cerr << "* Bin noch in [lin_dep]\n" << flush; |
---|
| 103 | |
---|
| 104 | // === Quest for a pivot === |
---|
| 105 | for ( i=1; i<=max_columns; i++) |
---|
[806c18] | 106 | if( A(now_row,i) != 0 ) |
---|
[6ead9d] | 107 | { |
---|
| 108 | pivot[now_row] = i; |
---|
| 109 | break;// i = max_columns; |
---|
[806c18] | 110 | |
---|
[6ead9d] | 111 | } |
---|
| 112 | // cout << "* pivot["<<now_row<<"] = " << pivot[now_row] << endl << flush; |
---|
| 113 | |
---|
| 114 | // === Is the reduced row (0,...,0)? ==== |
---|
| 115 | if( pivot[now_row] == 0 ) |
---|
| 116 | { |
---|
| 117 | if( INVERSE ) |
---|
| 118 | for ( i=1; i<=max_columns; i++) |
---|
[806c18] | 119 | A(now_row, max_columns+i) = 0; |
---|
[6ead9d] | 120 | return (true); |
---|
| 121 | } |
---|
| 122 | // === The row is linear independent ==== |
---|
| 123 | // cout << "* reduzierte Matrix (1) ist " << A << endl; |
---|
| 124 | if( INVERSE ) // Identity matrix |
---|
| 125 | A(now_row, now_row+max_columns) = 1; |
---|
| 126 | // cout << "* reduzierte Matrix (2) ist " << A << endl; |
---|
| 127 | return(false); |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | |
---|
| 131 | //--<>--------------------------------- |
---|
| 132 | void LGS::reduce(int fix, int row) |
---|
| 133 | //--<>--------------------------------- |
---|
| 134 | { // Reduce the `rowŽ by means of row `fixŽ |
---|
| 135 | |
---|
| 136 | if( A(row, pivot[fix]) == 0 ) return; |
---|
| 137 | |
---|
| 138 | CanonicalForm mul = A(row, pivot[fix]) / A(fix, pivot[fix]); |
---|
| 139 | |
---|
| 140 | // cout << "* Multiplikationsfaktor ist " << mul << endl; |
---|
[806c18] | 141 | // cout << "* Aktuelle Matrix ist " << A |
---|
[6ead9d] | 142 | // << "\n b = " << b << endl; |
---|
| 143 | // cout << "max_columns = "<< max_columns << endl; |
---|
| 144 | // cout << "* now_row = " << now_row <<", row=" << row << endl << flush; |
---|
| 145 | // cout << "* Pivot[" << row << "] = " << pivot[row] << endl << flush; |
---|
| 146 | // cout << "* Pivot[" << fix << "] = " << pivot[fix] << endl << flush; |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | for (int i=1; i<=max_columns; i++) |
---|
[806c18] | 150 | if( A(fix,i) != 0 ) |
---|
[6ead9d] | 151 | A(row, i) -= mul* A(fix,i); |
---|
[806c18] | 152 | if( b(fix,1) != 0 ) |
---|
[6ead9d] | 153 | b(row,1) -= mul * b(fix,1); |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | if ( INVERSE ) |
---|
| 157 | for (int i=max_columns+1; i<=2*max_columns; i++) |
---|
[806c18] | 158 | if( A(fix,i) != 0 ) |
---|
| 159 | A(row, i) -= (mul* A(fix,i)); |
---|
[6ead9d] | 160 | // A(row,i) = 777; |
---|
| 161 | |
---|
| 162 | // cout << "* reduzierte Matrix ist " << A << endl; |
---|
| 163 | // cout << "* Verlasse [reduce]\n" << flush; |
---|
| 164 | |
---|
| 165 | } |
---|
| 166 | /////////////////////////////////////////////////////// |
---|
| 167 | // The Matrix (A | I) has maximal rank and |
---|
| 168 | // the ŽAŽ part has to be normalized only |
---|
| 169 | //--<>--------------------------------- |
---|
| 170 | void LGS::inverse( CFMatrix & I ) |
---|
| 171 | //--<>--------------------------------- |
---|
[806c18] | 172 | { |
---|
[6ead9d] | 173 | int i,j; |
---|
| 174 | CanonicalForm mul; |
---|
| 175 | |
---|
| 176 | // cout << "* PRÄ-Invers:" << A << endl; |
---|
| 177 | for(i=1; i<=max_rows; i++) |
---|
| 178 | { |
---|
| 179 | for( j=max_columns+1; j<=2*max_columns; j++) |
---|
| 180 | if( A(i,j)!=0 ) |
---|
[806c18] | 181 | A(i,j) /= A( i, pivot[i]); |
---|
[6ead9d] | 182 | A( i, pivot[i]) = 1; |
---|
| 183 | } |
---|
| 184 | // cout << "* Inverse Matrix ist " << A << endl; |
---|
| 185 | |
---|
| 186 | for(i=1; i<=max_rows; i++) |
---|
| 187 | { |
---|
| 188 | // cout << "pivot["<<i<<"] = " << pivot[i]<< endl; |
---|
| 189 | for(j=1; j<=max_columns; j++) |
---|
| 190 | I(pivot[i],j) = A(i, j+max_columns); |
---|
| 191 | } |
---|
| 192 | } |
---|
| 193 | //--<>--------------------------------- |
---|
| 194 | int LGS::rank(void) |
---|
| 195 | //--<>--------------------------------- |
---|
| 196 | { // Return the current rank of the matrix |
---|
| 197 | return( now_row-1); |
---|
| 198 | } |
---|
| 199 | |
---|
| 200 | |
---|
| 201 | //--<>--------------------------------- |
---|
| 202 | void LGS::print(void) |
---|
| 203 | //--<>--------------------------------- |
---|
| 204 | { // Return the current rank of the matrix |
---|
[21d54bf] | 205 | #ifndef NOSTREAMIO |
---|
[6ead9d] | 206 | cout << "A = " << A << "\nb = " << b << endl; |
---|
| 207 | #endif |
---|
| 208 | } |
---|
| 209 | |
---|
| 210 | |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | //--<>--------------------------------- |
---|
| 214 | int LGS::corank(void) |
---|
| 215 | //--<>--------------------------------- |
---|
| 216 | { // Return the current rank of the matrix |
---|
| 217 | return( ( (max_rows < max_columns ) ? max_rows : max_columns )-now_row+1 ); |
---|
| 218 | } |
---|
| 219 | |
---|
| 220 | //--<>--------------------------------- |
---|
| 221 | int LGS::ErgCol(int row, int basis[]) |
---|
| 222 | //--<>--------------------------------- |
---|
[806c18] | 223 | { |
---|
[6ead9d] | 224 | bool state = false; |
---|
| 225 | |
---|
| 226 | for( int i=1; i<=max_columns; i++) |
---|
| 227 | { |
---|
| 228 | if( A(row,i) != 0 ) |
---|
| 229 | { |
---|
| 230 | state = true; |
---|
| 231 | for( int j=1; j<=basis[0]; j++) |
---|
[806c18] | 232 | if( i==basis[j] ) |
---|
| 233 | state = false; |
---|
| 234 | if( state == true ) |
---|
| 235 | return (i); |
---|
[6ead9d] | 236 | } |
---|
| 237 | } |
---|
| 238 | // A row contains only pivot entry/-ies |
---|
| 239 | for( int j=1; j<=basis[0]; j++) |
---|
| 240 | if( A(row, basis[j]) != 0 ) |
---|
| 241 | return( -basis[j] ); |
---|
| 242 | |
---|
[806c18] | 243 | |
---|
[6ead9d] | 244 | AUSGABE_LGS("Zeile ist " << row << endl); |
---|
| 245 | AUSGABE_ERR("* Mistake in [lgs.cc]! Impossible result. Aborting!\n"); |
---|
| 246 | exit (1); |
---|
| 247 | } |
---|
| 248 | |
---|
| 249 | |
---|
| 250 | //--<>--------------------------------- |
---|
| 251 | CFMatrix LGS::GetKernelBasis(void) |
---|
| 252 | //--<>--------------------------------- |
---|
| 253 | { |
---|
| 254 | int i,z; |
---|
[806c18] | 255 | int dim = corank(); |
---|
[6ead9d] | 256 | |
---|
| 257 | bool* tmp_vec = new bool[max_columns+1]; |
---|
| 258 | int* basis = new int [dim+1]; |
---|
| 259 | basis[0]=1; |
---|
| 260 | CFMatrix erg (dim, max_columns); // Each row is one solution |
---|
[806c18] | 261 | |
---|
[6ead9d] | 262 | // === Searching free parameters === |
---|
[806c18] | 263 | for( i=1; i<=max_rows; i++) |
---|
[6ead9d] | 264 | tmp_vec[i]=false; |
---|
| 265 | for( i=1; i<now_row; i++) |
---|
[806c18] | 266 | if( pivot[i] != 0 ) |
---|
[6ead9d] | 267 | tmp_vec[pivot[i]] = true; |
---|
| 268 | |
---|
| 269 | for( i=1; i<=max_columns; i++) |
---|
| 270 | if(tmp_vec[i] == false) |
---|
| 271 | { |
---|
| 272 | basis[basis[0]]=i; |
---|
| 273 | basis[0]++; |
---|
| 274 | } |
---|
| 275 | ASSERT(basis[0]-1==dim, "Wrong dimensions"); |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | // cout << "* Freie Parameter sind:"; |
---|
| 279 | // for ( int k=1; k<basis[0]; k++) cout << basis[k] << ", "; |
---|
| 280 | // cout << endl; |
---|
| 281 | |
---|
| 282 | // === Creating the basis vectors === |
---|
| 283 | for(z=1; z<now_row; z++) |
---|
| 284 | for(i=1; i<=dim; i++) |
---|
| 285 | { |
---|
| 286 | erg(i, basis[i]) = 1; |
---|
| 287 | erg(i, pivot[z]) = -(A(z, basis[i]) / A(z,pivot[z])); |
---|
| 288 | } |
---|
[806c18] | 289 | |
---|
[6ead9d] | 290 | // cout << "Kernbasis ist " << erg; |
---|
| 291 | return erg; |
---|
| 292 | } |
---|
| 293 | |
---|
| 294 | //--<>--------------------------------- |
---|
| 295 | CFMatrix LGS::GetSolutionVector(void) |
---|
| 296 | //--<>--------------------------------- |
---|
| 297 | { // Ax=b has exactly one solution vector x |
---|
| 298 | int z; |
---|
[806c18] | 299 | CFMatrix erg (1, max_columns); |
---|
| 300 | |
---|
[6ead9d] | 301 | // === Creating the basis vectors === |
---|
| 302 | for(z=1; z<=max_columns; z++) |
---|
| 303 | { |
---|
| 304 | erg(1,pivot[z]) = b(z,1); // Vector b |
---|
| 305 | erg(1,pivot[z]) /= A(z,pivot[z]); |
---|
| 306 | } |
---|
| 307 | return erg; |
---|
| 308 | } |
---|
| 309 | #endif |
---|