[35aab3] | 1 | // ---------------------------------------------------------------------------- |
---|
| 2 | // kmatrix.h |
---|
| 3 | // begin of file |
---|
| 4 | // Stephan Endrass, endrass@mathematik.uni-mainz.de |
---|
| 5 | // 23.7.99 |
---|
| 6 | // ---------------------------------------------------------------------------- |
---|
| 7 | |
---|
| 8 | #ifndef KMATRIX_H |
---|
| 9 | #define KMATRIX_H |
---|
| 10 | |
---|
| 11 | // ---------------------------------------------------------------------------- |
---|
| 12 | // template class for matrices with coefficients in the field K |
---|
| 13 | // K is a class representing elements of a field |
---|
| 14 | // The implementation of K is expected to have overloaded |
---|
| 15 | // the operators +, -, *, /, +=, -=, *= and /=. |
---|
| 16 | // The expressions (K)0 and (K)1 should cast to the 0 and 1 of K. |
---|
| 17 | // Additionally we use the following functions in class K: |
---|
| 18 | // |
---|
| 19 | // member functions: |
---|
| 20 | // |
---|
| 21 | // double complexity( void ); |
---|
| 22 | // |
---|
| 23 | // friend functions: |
---|
| 24 | // |
---|
| 25 | // friend K gcd( const K &a,const K &b ); // gcd(a,b) |
---|
| 26 | // friend K gcd( K* a,int k ); // gcd(a[0],...,a[k-1]) |
---|
| 27 | // |
---|
| 28 | // The complexity function should return a measure indicating |
---|
| 29 | // how complicated this number is in terms of memory usage |
---|
| 30 | // and arithmetic operations. For a rational p/q, one could |
---|
| 31 | // return max(|p|,|q|). This fuction is used for pivoting. |
---|
| 32 | // |
---|
| 33 | // The gcd of two numbers a,b should be a number g such that |
---|
| 34 | // the complexities of a/g and b/g are less or equal than those |
---|
| 35 | // of a and b. For rationals p1/q1, p2/q2 one could return the |
---|
| 36 | // quotient of integer gcd's gcd(p1,p2)/gcd(q1,q2). |
---|
| 37 | // |
---|
| 38 | // ---------------------------------------------------------------------------- |
---|
| 39 | |
---|
| 40 | template<class K> class KMatrix |
---|
| 41 | { |
---|
| 42 | |
---|
| 43 | private: |
---|
| 44 | |
---|
| 45 | K *a; // the entries ot the matrix |
---|
| 46 | int rows; // number of rows |
---|
| 47 | int cols; // number of columns |
---|
| 48 | |
---|
| 49 | public: |
---|
| 50 | |
---|
| 51 | KMatrix( ); // init zero |
---|
| 52 | KMatrix( const KMatrix& ); // copy constructor |
---|
| 53 | KMatrix( int,int ); // preallocate rows & columns |
---|
| 54 | |
---|
| 55 | ~KMatrix( ); // destructor |
---|
| 56 | |
---|
| 57 | void copy_delete ( void ); // delete associated memory |
---|
| 58 | void copy_new ( int ); // allocate associated memory |
---|
| 59 | void copy_zero ( void ); // init zero |
---|
| 60 | void copy_unit ( int ); // init as unit matrix |
---|
| 61 | void copy_shallow( KMatrix& ); // shallow copy |
---|
| 62 | void copy_deep ( const KMatrix& ); // deep copy |
---|
| 63 | |
---|
| 64 | K get( int,int ) const; // get an element |
---|
| 65 | void set( int,int,const K& ); // set an element |
---|
| 66 | |
---|
| 67 | int row_is_zero( int ) const; // test if row is zero |
---|
| 68 | int column_is_zero( int ) const; // test if column is zero |
---|
| 69 | |
---|
| 70 | int column_pivot( int,int ) const; |
---|
| 71 | |
---|
| 72 | int gausseliminate( void ); // Gauss elimination |
---|
| 73 | int rank( void ) const; // compute the rank |
---|
| 74 | int solve( K**,int* ); // solve Ax=b from (A|b) |
---|
| 75 | |
---|
| 76 | // elementary transformations |
---|
| 77 | |
---|
| 78 | K multiply_row( int,const K& ); |
---|
| 79 | K add_rows( int,int,const K&,const K& ); |
---|
| 80 | int swap_rows( int,int ); |
---|
| 81 | K set_row_primitive( int ); |
---|
| 82 | |
---|
| 83 | int is_quadratic( void ) const; |
---|
| 84 | int is_symmetric( void ) const; |
---|
| 85 | |
---|
| 86 | K determinant( void ) const; |
---|
| 87 | |
---|
| 88 | #ifdef KMATRIX_DEBUG |
---|
| 89 | void test_row( int ) const; |
---|
| 90 | void test_col( int ) const; |
---|
| 91 | #endif |
---|
| 92 | |
---|
| 93 | #ifdef KMATRIX_PRINT |
---|
| 94 | friend ostream & operator << ( ostream&,const KMatrix& ); |
---|
| 95 | #endif |
---|
| 96 | }; |
---|
| 97 | |
---|
| 98 | // ------------------------------------ |
---|
| 99 | // inline functions for class KMatrix |
---|
| 100 | // ------------------------------------ |
---|
| 101 | |
---|
| 102 | // ---------------------------------------------------------------------------- |
---|
| 103 | // Delete memory associated to a KMatrix |
---|
| 104 | // ---------------------------------------------------------------------------- |
---|
| 105 | |
---|
| 106 | template<class K> |
---|
| 107 | inline void KMatrix<K>::copy_delete( void ) |
---|
| 108 | { |
---|
| 109 | if( a != (K*)NULL && rows > 0 && cols > 0 ) delete [] a; |
---|
| 110 | copy_zero( ); |
---|
| 111 | } |
---|
| 112 | |
---|
| 113 | // ---------------------------------------------------------------------------- |
---|
| 114 | // Allocate memory associated to a KMatrix |
---|
| 115 | // ---------------------------------------------------------------------------- |
---|
| 116 | |
---|
| 117 | template<class K> |
---|
| 118 | inline void KMatrix<K>::copy_new( int k ) |
---|
| 119 | { |
---|
| 120 | if( k > 0 ) |
---|
| 121 | { |
---|
| 122 | a = new K[k]; |
---|
| 123 | |
---|
| 124 | #ifndef NDEBUG |
---|
| 125 | if( a == (K*)NULL ) |
---|
| 126 | { |
---|
| 127 | #ifdef KMATRIX_PRINT |
---|
| 128 | #ifdef KMATRIX_IOSTREAM |
---|
| 129 | cerr << "void KMatrix::copy_new( int k )"; |
---|
| 130 | cerr << ": no memory left ..." << endl; |
---|
| 131 | #else |
---|
| 132 | fprintf( stderr,"void KMatrix::copy_new( int k )" ); |
---|
| 133 | fprintf( stderr,": no memory left ...\n" ); |
---|
| 134 | #endif |
---|
| 135 | #endif |
---|
| 136 | |
---|
| 137 | exit( 1 ); |
---|
| 138 | } |
---|
| 139 | #endif |
---|
| 140 | } |
---|
| 141 | else if( k == 0 ) |
---|
| 142 | { |
---|
| 143 | a = (K*)NULL; |
---|
| 144 | } |
---|
| 145 | else |
---|
| 146 | { |
---|
| 147 | #ifdef KMATRIX_PRINT |
---|
| 148 | #ifdef KMATRIX_IOSTREAM |
---|
| 149 | cerr << "void KMatrix::copy_new( int k )"; |
---|
| 150 | cerr << ": k < 0 ..." << endl; |
---|
| 151 | #else |
---|
| 152 | fprintf( stderr,"void KMatrix::copy_new( int k )" ); |
---|
| 153 | fprintf( stderr,": k < 0 ...\n" ); |
---|
| 154 | #endif |
---|
| 155 | #endif |
---|
| 156 | |
---|
| 157 | exit( 1 ); |
---|
| 158 | } |
---|
| 159 | } |
---|
| 160 | |
---|
| 161 | // ---------------------------------------------------------------------------- |
---|
| 162 | // Initialize a KMatrix with 0 |
---|
| 163 | // ---------------------------------------------------------------------------- |
---|
| 164 | |
---|
| 165 | template<class K> |
---|
| 166 | inline void KMatrix<K>::copy_zero( void ) |
---|
| 167 | { |
---|
| 168 | a = (K*)NULL; |
---|
| 169 | rows = cols = 0; |
---|
| 170 | } |
---|
| 171 | |
---|
| 172 | // ---------------------------------------------------------------------------- |
---|
| 173 | // Initialize a KMatrix with the unit matrix |
---|
| 174 | // ---------------------------------------------------------------------------- |
---|
| 175 | |
---|
| 176 | template<class K> |
---|
| 177 | inline void KMatrix<K>::copy_unit( int rank ) |
---|
| 178 | { |
---|
| 179 | int r,n=rank*rank; |
---|
| 180 | copy_new( n ); |
---|
| 181 | rows = cols = rank; |
---|
| 182 | |
---|
| 183 | for( r=0; r<n; a[r++]=(K)0 ); |
---|
| 184 | |
---|
| 185 | for( r=0; r<rows; r++ ) |
---|
| 186 | { |
---|
| 187 | a[r*cols+r] = (K)1; |
---|
| 188 | } |
---|
| 189 | } |
---|
| 190 | |
---|
| 191 | // ---------------------------------------------------------------------------- |
---|
| 192 | // Shallow copy |
---|
| 193 | // ---------------------------------------------------------------------------- |
---|
| 194 | |
---|
| 195 | template<class K> |
---|
| 196 | inline void KMatrix<K>::copy_shallow( KMatrix &m ) |
---|
| 197 | { |
---|
| 198 | a = m.a; |
---|
| 199 | rows = m.rows; |
---|
| 200 | cols = m.cols; |
---|
| 201 | } |
---|
| 202 | |
---|
| 203 | // ---------------------------------------------------------------------------- |
---|
| 204 | // Deep copy |
---|
| 205 | // ---------------------------------------------------------------------------- |
---|
| 206 | |
---|
| 207 | template<class K> |
---|
| 208 | inline void KMatrix<K>::copy_deep( const KMatrix &m ) |
---|
| 209 | { |
---|
| 210 | if( m.a == (K*)NULL ) |
---|
| 211 | { |
---|
| 212 | copy_zero( ); |
---|
| 213 | } |
---|
| 214 | else |
---|
| 215 | { |
---|
| 216 | int n=m.rows*m.cols; |
---|
| 217 | copy_new( n ); |
---|
| 218 | rows = m.rows; |
---|
| 219 | cols = m.cols; |
---|
| 220 | |
---|
| 221 | for( int i=0; i<n; i++ ) |
---|
| 222 | { |
---|
| 223 | a[i] = m.a[i]; |
---|
| 224 | } |
---|
| 225 | } |
---|
| 226 | } |
---|
| 227 | |
---|
| 228 | // ---------------------------------------------------------------------------- |
---|
| 229 | // Zero constructor |
---|
| 230 | // ---------------------------------------------------------------------------- |
---|
| 231 | |
---|
| 232 | template<class K> |
---|
| 233 | inline KMatrix<K>::KMatrix( ) |
---|
| 234 | { |
---|
| 235 | copy_zero( ); |
---|
| 236 | } |
---|
| 237 | |
---|
| 238 | // ---------------------------------------------------------------------------- |
---|
| 239 | // Copy constructor |
---|
| 240 | // ---------------------------------------------------------------------------- |
---|
| 241 | |
---|
| 242 | template<class K> |
---|
| 243 | inline KMatrix<K>::KMatrix( const KMatrix &m ) |
---|
| 244 | { |
---|
| 245 | copy_deep( m ); |
---|
| 246 | } |
---|
| 247 | |
---|
| 248 | // ---------------------------------------------------------------------------- |
---|
| 249 | // Zero r by c matrix constructor |
---|
| 250 | // ---------------------------------------------------------------------------- |
---|
| 251 | |
---|
| 252 | template<class K> |
---|
| 253 | KMatrix<K>::KMatrix( int r,int c ) |
---|
| 254 | { |
---|
| 255 | int n = r*c; |
---|
| 256 | |
---|
| 257 | copy_new( n ); |
---|
| 258 | |
---|
| 259 | rows = r; |
---|
| 260 | cols = c; |
---|
| 261 | |
---|
| 262 | for( int i=0; i<n; i++ ) |
---|
| 263 | { |
---|
| 264 | a[i]=(K)0; |
---|
| 265 | } |
---|
| 266 | } |
---|
| 267 | |
---|
| 268 | // ---------------------------------------------------------------------------- |
---|
| 269 | // Destructor |
---|
| 270 | // ---------------------------------------------------------------------------- |
---|
| 271 | |
---|
| 272 | template<class K> |
---|
| 273 | KMatrix<K>::~KMatrix( ) |
---|
| 274 | { |
---|
| 275 | copy_delete( ); |
---|
| 276 | } |
---|
| 277 | |
---|
| 278 | // ------------------------------------------------- |
---|
| 279 | // non-inline template functions for class KMatrix |
---|
| 280 | // ------------------------------------------------- |
---|
| 281 | |
---|
| 282 | // ---------------------------------------------------------------------------- |
---|
| 283 | // Debugging functions |
---|
| 284 | // ---------------------------------------------------------------------------- |
---|
| 285 | |
---|
| 286 | #ifdef KMATRIX_DEBUG |
---|
| 287 | |
---|
| 288 | template<class K> |
---|
| 289 | void KMatrix<K>::test_row( int r ) const |
---|
| 290 | { |
---|
| 291 | if( r<0 || r>=rows ) |
---|
| 292 | { |
---|
| 293 | #ifdef KMATRIX_PRINT |
---|
| 294 | #ifdef KMATRIX_IOSTREAM |
---|
| 295 | cerr << "KMatrix<K>::test_row( " << r << " )" << endl; |
---|
| 296 | cerr << " rows = " << rows << endl; |
---|
| 297 | cerr << " exiting...." << endl; |
---|
| 298 | #else |
---|
| 299 | fprintf( stderr,"KMatrix<K>::test_row( %d )\n",r ); |
---|
| 300 | fprintf( stderr," rows = %d\n",rows ); |
---|
| 301 | fprintf( stderr," exiting....\n" ); |
---|
| 302 | #endif |
---|
| 303 | #endif |
---|
| 304 | |
---|
| 305 | exit( 1 ); |
---|
| 306 | } |
---|
| 307 | } |
---|
| 308 | |
---|
| 309 | template<class K> |
---|
| 310 | void KMatrix<K>::test_col( int c ) const |
---|
| 311 | { |
---|
| 312 | if( c<0 || c>=cols ) |
---|
| 313 | { |
---|
| 314 | #ifdef KMATRIX_PRINT |
---|
| 315 | #ifdef KMATRIX_IOSTREAM |
---|
| 316 | cerr << "KMatrix<K>::test_col( " << c << " )" << endl; |
---|
| 317 | cerr << " cols = " << cols << endl; |
---|
| 318 | cerr << " exiting...." << endl; |
---|
| 319 | #else |
---|
| 320 | fprintf( stderr,"KMatrix<K>::test_col( %d )\n",c ); |
---|
| 321 | fprintf( stderr," cols = %d\n",cols ); |
---|
| 322 | fprintf( stderr," exiting....\n" ); |
---|
| 323 | #endif |
---|
| 324 | #endif |
---|
| 325 | |
---|
| 326 | exit( 1 ); |
---|
| 327 | } |
---|
| 328 | } |
---|
| 329 | |
---|
| 330 | #endif |
---|
| 331 | |
---|
| 332 | // ---------------------------------------------------------------------------- |
---|
| 333 | // get coefficient at row r and column c |
---|
| 334 | // return value: the coefficient |
---|
| 335 | // ---------------------------------------------------------------------------- |
---|
| 336 | |
---|
| 337 | template<class K> |
---|
| 338 | K KMatrix<K>::get( int r,int c ) const |
---|
| 339 | { |
---|
| 340 | #ifdef KMATRIX_DEBUG |
---|
| 341 | test_row( r ); |
---|
| 342 | test_col( c ); |
---|
| 343 | #endif |
---|
| 344 | |
---|
| 345 | return a[r*cols+c]; |
---|
| 346 | } |
---|
| 347 | |
---|
| 348 | // ---------------------------------------------------------------------------- |
---|
| 349 | // sets coefficient at row r and column c to value |
---|
| 350 | // ---------------------------------------------------------------------------- |
---|
| 351 | |
---|
| 352 | template<class K> |
---|
| 353 | void KMatrix<K>::set( int r,int c,const K &value ) |
---|
| 354 | { |
---|
| 355 | #ifdef KMATRIX_DEBUG |
---|
| 356 | test_row( r ); |
---|
| 357 | test_col( c ); |
---|
| 358 | #endif |
---|
| 359 | |
---|
| 360 | a[r*cols+c] = value; |
---|
| 361 | } |
---|
| 362 | |
---|
| 363 | // ---------------------------------------------------------------------------- |
---|
| 364 | // interchanges the rows r1 and r2 |
---|
| 365 | // return value: 1 if r1==r2 |
---|
| 366 | // return value: -1 if r1!=r2 |
---|
| 367 | // caution: the determinant changes its sign by the return value |
---|
| 368 | // ---------------------------------------------------------------------------- |
---|
| 369 | |
---|
| 370 | template<class K> |
---|
| 371 | int KMatrix<K>::swap_rows( int r1,int r2 ) |
---|
| 372 | { |
---|
| 373 | #ifdef KMATRIX_DEBUG |
---|
| 374 | test_row( r1 ); |
---|
| 375 | test_row( r2 ); |
---|
| 376 | #endif |
---|
| 377 | |
---|
| 378 | if( r1 == r2 ) return 1; |
---|
| 379 | |
---|
| 380 | K tmp; |
---|
| 381 | |
---|
| 382 | for( int c=0; c<cols; c++ ) |
---|
| 383 | { |
---|
| 384 | tmp = a[r1*cols+c]; |
---|
| 385 | a[r1*cols+c] = a[r2*cols+c]; |
---|
| 386 | a[r2*cols+c] = tmp; |
---|
| 387 | } |
---|
| 388 | |
---|
| 389 | return -1; |
---|
| 390 | } |
---|
| 391 | |
---|
| 392 | // ---------------------------------------------------------------------------- |
---|
| 393 | // replaces row r by its multiple (row r)*factor |
---|
| 394 | // return value: factor |
---|
| 395 | // caution: the determinant changes by the return value |
---|
| 396 | // ---------------------------------------------------------------------------- |
---|
| 397 | |
---|
| 398 | template<class K> |
---|
| 399 | K KMatrix<K>::multiply_row( int r,const K &factor ) |
---|
| 400 | { |
---|
| 401 | #ifdef KMATRIX_DEBUG |
---|
| 402 | test_row( r ); |
---|
| 403 | #endif |
---|
| 404 | |
---|
| 405 | int i_src = r*cols; |
---|
| 406 | |
---|
| 407 | for( int i=0; i<cols; i++,i_src++ ) |
---|
| 408 | { |
---|
| 409 | a[i_src] *= factor; |
---|
| 410 | } |
---|
| 411 | return factor; |
---|
| 412 | } |
---|
| 413 | |
---|
| 414 | // ---------------------------------------------------------------------------- |
---|
| 415 | // replaces row dest by the linear combination |
---|
| 416 | // (row src)*factor_src + (row dest)*factor_dest |
---|
| 417 | // return value: factor_dest |
---|
| 418 | // caution: the determinant changes by the return value |
---|
| 419 | // ---------------------------------------------------------------------------- |
---|
| 420 | |
---|
| 421 | template<class K> |
---|
| 422 | K KMatrix<K>::add_rows( |
---|
| 423 | int src,int dest,const K &factor_src,const K &factor_dest ) |
---|
| 424 | { |
---|
| 425 | #ifdef KMATRIX_DEBUG |
---|
| 426 | test_row( src ); |
---|
| 427 | test_row( dest ); |
---|
| 428 | #endif |
---|
| 429 | |
---|
| 430 | int i; |
---|
| 431 | int i_src = src*cols; |
---|
| 432 | int i_dest = dest*cols; |
---|
| 433 | |
---|
| 434 | for( i=0; i<cols; i++,i_src++,i_dest++ ) |
---|
| 435 | { |
---|
| 436 | a[i_dest] = a[i_src]*factor_src + a[i_dest]*factor_dest; |
---|
| 437 | } |
---|
| 438 | |
---|
| 439 | return factor_dest; |
---|
| 440 | } |
---|
| 441 | |
---|
| 442 | // ---------------------------------------------------------------------------- |
---|
| 443 | // test if row r is zero |
---|
| 444 | // return value: TRUE if zero |
---|
| 445 | // FALSE if not zero |
---|
| 446 | // ---------------------------------------------------------------------------- |
---|
| 447 | |
---|
| 448 | template<class K> |
---|
| 449 | int KMatrix<K>::row_is_zero( int r ) const |
---|
| 450 | { |
---|
| 451 | #ifdef KMATRIX_DEBUG |
---|
| 452 | test_row( r ); |
---|
| 453 | #endif |
---|
| 454 | |
---|
| 455 | for( int c=0; c<cols; c++ ) |
---|
| 456 | { |
---|
| 457 | if( a[r*cols+c] != (K)0 ) return FALSE; |
---|
| 458 | } |
---|
| 459 | return TRUE; |
---|
| 460 | } |
---|
| 461 | |
---|
| 462 | // ---------------------------------------------------------------------------- |
---|
| 463 | // test if column c is zero |
---|
| 464 | // return value: TRUE if zero |
---|
| 465 | // FALSE if not zero |
---|
| 466 | // ---------------------------------------------------------------------------- |
---|
| 467 | |
---|
| 468 | template<class K> |
---|
| 469 | int KMatrix<K>::column_is_zero( int c ) const |
---|
| 470 | { |
---|
| 471 | #ifdef KMATRIX_DEBUG |
---|
| 472 | test_col( c ); |
---|
| 473 | #endif |
---|
| 474 | |
---|
| 475 | for( int r=0; r<rows; r++ ) |
---|
| 476 | { |
---|
| 477 | if( a[r*cols+c] != (K)0 ) return FALSE; |
---|
| 478 | } |
---|
| 479 | return TRUE; |
---|
| 480 | } |
---|
| 481 | |
---|
| 482 | // ---------------------------------------------------------------------------- |
---|
| 483 | // find the element of column c if smallest nonzero absolute value |
---|
| 484 | // consider only elements in row r0 or below |
---|
| 485 | // return value: the row of the element |
---|
| 486 | // ---------------------------------------------------------------------------- |
---|
| 487 | |
---|
| 488 | template<class K> |
---|
| 489 | int KMatrix<K>::column_pivot( int r0,int c ) const |
---|
| 490 | { |
---|
| 491 | #ifdef KMATRIX_DEBUG |
---|
| 492 | test_row( r0 ); |
---|
| 493 | test_col( c ); |
---|
| 494 | #endif |
---|
| 495 | |
---|
| 496 | int r; |
---|
| 497 | // find first nonzero entry in column c |
---|
| 498 | |
---|
| 499 | for( r=r0; r<rows && a[r*cols+c]==(K)0; r++ ); |
---|
| 500 | |
---|
| 501 | if( r == rows ) |
---|
| 502 | { |
---|
| 503 | // column is zero |
---|
| 504 | |
---|
| 505 | return -1; |
---|
| 506 | } |
---|
| 507 | else |
---|
| 508 | { |
---|
| 509 | double val = a[r*cols+c].complexity( ); |
---|
| 510 | double val_new = 0.0; |
---|
| 511 | int pivot = r; |
---|
| 512 | |
---|
| 513 | for( ; r<rows; r++ ) |
---|
| 514 | { |
---|
| 515 | if( a[r*cols+c] != (K)0 && |
---|
| 516 | ( val_new = a[r*cols+c].complexity( ) ) < val ) |
---|
| 517 | { |
---|
| 518 | val = val_new; |
---|
| 519 | pivot = r; |
---|
| 520 | } |
---|
| 521 | } |
---|
| 522 | return pivot; |
---|
| 523 | } |
---|
| 524 | } |
---|
| 525 | |
---|
| 526 | // ---------------------------------------------------------------------------- |
---|
| 527 | // divide row r by the gcd of all elements |
---|
| 528 | // ---------------------------------------------------------------------------- |
---|
| 529 | |
---|
| 530 | template<class K> |
---|
| 531 | K KMatrix<K>::set_row_primitive( int r ) |
---|
| 532 | { |
---|
| 533 | #ifdef KMATRIX_DEBUG |
---|
| 534 | test_row( r ); |
---|
| 535 | #endif |
---|
| 536 | |
---|
| 537 | K g = gcd( &(a[r*cols]),cols ); |
---|
| 538 | |
---|
| 539 | for( int c=0; c<cols; c++ ) |
---|
| 540 | { |
---|
| 541 | a[r*cols+c] /= g; |
---|
| 542 | } |
---|
| 543 | return g; |
---|
| 544 | } |
---|
| 545 | |
---|
| 546 | // ---------------------------------------------------------------------------- |
---|
| 547 | // convert the matrix to upper triangular form |
---|
| 548 | // return value: rank of the matrix |
---|
| 549 | // ---------------------------------------------------------------------------- |
---|
| 550 | |
---|
| 551 | template<class K> |
---|
| 552 | int KMatrix<K>::gausseliminate( void ) |
---|
| 553 | { |
---|
| 554 | int r,c,rank = 0; |
---|
| 555 | K g; |
---|
| 556 | |
---|
| 557 | // make sure that the elements of each row have gcd=1 |
---|
| 558 | // this is useful for pivoting |
---|
| 559 | |
---|
| 560 | for( r=0; r<rows; r++ ) |
---|
| 561 | { |
---|
| 562 | set_row_primitive( r ); |
---|
| 563 | } |
---|
| 564 | |
---|
| 565 | // search a pivoting element in each column |
---|
| 566 | // perform Gauss elimination |
---|
| 567 | |
---|
| 568 | for( c=0; c<cols && rank<rows; c++ ) |
---|
| 569 | { |
---|
| 570 | if( ( r = column_pivot( rank,c )) >= 0 ) |
---|
| 571 | { |
---|
| 572 | swap_rows( rank,r ); |
---|
| 573 | |
---|
| 574 | for( r=rank+1; r<rows; r++ ) |
---|
| 575 | { |
---|
| 576 | if( a[r*cols+c] != (K)0 ) |
---|
| 577 | { |
---|
| 578 | g = gcd( a[r*cols+c],a[rank*cols+c] ); |
---|
| 579 | add_rows( rank,r,-a[r*cols+c]/g,a[rank*cols+c]/g ); |
---|
| 580 | set_row_primitive( r ); |
---|
| 581 | } |
---|
| 582 | } |
---|
| 583 | rank++; |
---|
| 584 | } |
---|
| 585 | } |
---|
| 586 | return rank; |
---|
| 587 | } |
---|
| 588 | |
---|
| 589 | // ---------------------------------------------------------------------------- |
---|
| 590 | // solve the linear system of equations given by |
---|
| 591 | // (x1,...,xn,-1)*(*this) = 0 |
---|
| 592 | // return value: rank of the matrix |
---|
| 593 | // k is set to the number of variables |
---|
| 594 | // rat[0],...,rat[k-1] are set to the solutions |
---|
| 595 | // ---------------------------------------------------------------------------- |
---|
| 596 | |
---|
| 597 | template<class K> |
---|
| 598 | int KMatrix<K>::solve( K **solution,int *k ) |
---|
| 599 | { |
---|
| 600 | int r,c,rank = 0; |
---|
| 601 | K g; |
---|
| 602 | |
---|
| 603 | // ---------------------------------------------------- |
---|
| 604 | // make sure that the elements of each row have gcd=1 |
---|
| 605 | // this is useful for pivoting |
---|
| 606 | // ---------------------------------------------------- |
---|
| 607 | |
---|
| 608 | for( r=0; r<rows; r++ ) |
---|
| 609 | { |
---|
| 610 | set_row_primitive( r ); |
---|
| 611 | } |
---|
| 612 | |
---|
| 613 | // ------------------------------------------ |
---|
| 614 | // search a pivoting element in each column |
---|
| 615 | // perform Gauss elimination |
---|
| 616 | // ------------------------------------------ |
---|
| 617 | |
---|
| 618 | for( c=0; c<cols && rank < rows; c++ ) |
---|
| 619 | { |
---|
| 620 | if( ( r = column_pivot( rank,c )) >= 0 ) |
---|
| 621 | { |
---|
| 622 | swap_rows( rank,r ); |
---|
| 623 | |
---|
| 624 | for( r=0; r<rank; r++ ) |
---|
| 625 | { |
---|
| 626 | if( a[r*cols+c] != (K)0 ) |
---|
| 627 | { |
---|
| 628 | g = gcd( a[r*cols+c],a[rank*cols+c] ); |
---|
| 629 | add_rows( rank,r,-a[r*cols+c]/g,a[rank*cols+c]/g ); |
---|
| 630 | set_row_primitive( r ); |
---|
| 631 | } |
---|
| 632 | } |
---|
| 633 | |
---|
| 634 | for( r=rank+1; r<rows; r++ ) |
---|
| 635 | { |
---|
| 636 | if( a[r*cols+c] != (K)0 ) |
---|
| 637 | { |
---|
| 638 | g = gcd( a[r*cols+c],a[rank*cols+c] ); |
---|
| 639 | add_rows( rank,r,-a[r*cols+c]/g,a[rank*cols+c]/g ); |
---|
| 640 | set_row_primitive( r ); |
---|
| 641 | } |
---|
| 642 | } |
---|
| 643 | |
---|
| 644 | rank++; |
---|
| 645 | } |
---|
| 646 | } |
---|
| 647 | |
---|
| 648 | if( rank < cols ) |
---|
| 649 | { |
---|
| 650 | // ---------------------- |
---|
| 651 | // equation is solvable |
---|
| 652 | // copy solutions |
---|
| 653 | // ---------------------- |
---|
| 654 | |
---|
| 655 | *solution = new K[cols-1]; |
---|
| 656 | *k = cols - 1; |
---|
| 657 | |
---|
| 658 | for( c=0; c<cols-1; c++ ) |
---|
| 659 | { |
---|
| 660 | (*solution)[c] = (K)0; |
---|
| 661 | } |
---|
| 662 | |
---|
| 663 | for( r=0; r<rows; r++ ) |
---|
| 664 | { |
---|
| 665 | for( c=0; c<cols && a[r*cols+c] == (K)0; c++ ); |
---|
| 666 | |
---|
| 667 | if( c < cols-1 ) |
---|
| 668 | { |
---|
| 669 | (*solution)[c] = ((K)a[(r+1)*cols-1])/a[r*cols+c]; |
---|
| 670 | } |
---|
| 671 | } |
---|
| 672 | } |
---|
| 673 | else |
---|
| 674 | { |
---|
| 675 | // -------------------------- |
---|
| 676 | // equation is not solvable |
---|
| 677 | // -------------------------- |
---|
| 678 | |
---|
| 679 | *solution = (K*)NULL; |
---|
| 680 | *k = 0; |
---|
| 681 | } |
---|
| 682 | |
---|
| 683 | return rank; |
---|
| 684 | } |
---|
| 685 | |
---|
| 686 | // ---------------------------------------------------------------------------- |
---|
| 687 | // compute the rank of the matrix |
---|
| 688 | // return value: rank of the matrix |
---|
| 689 | // ---------------------------------------------------------------------------- |
---|
| 690 | |
---|
| 691 | template<class K> |
---|
| 692 | int KMatrix<K>::rank( void ) const |
---|
| 693 | { |
---|
| 694 | KMatrix<K> dummy( *this ); |
---|
| 695 | |
---|
| 696 | return dummy.gausseliminate( ); |
---|
| 697 | } |
---|
| 698 | |
---|
| 699 | // ---------------------------------------------------------------------------- |
---|
| 700 | // print the matrix |
---|
| 701 | // return value: the output stream used |
---|
| 702 | // ---------------------------------------------------------------------------- |
---|
| 703 | |
---|
| 704 | #ifdef KMATRIX_PRINT |
---|
| 705 | |
---|
| 706 | template<class K> |
---|
| 707 | static |
---|
| 708 | void print_rational( ostream &s,int digits,const K &n ) |
---|
| 709 | { |
---|
| 710 | unsigned int num = digits - n.length( ); |
---|
| 711 | |
---|
| 712 | for( unsigned int i=0; i < num; i++ ) |
---|
| 713 | { |
---|
| 714 | #ifdef KMATRIX_IOSTREAM |
---|
| 715 | s << " "; |
---|
| 716 | #else |
---|
| 717 | fprintf( stdout," " ); |
---|
| 718 | #endif |
---|
| 719 | } |
---|
| 720 | |
---|
| 721 | s << n; |
---|
| 722 | } |
---|
| 723 | |
---|
| 724 | template<class K> |
---|
| 725 | ostream & operator << ( ostream &s,const KMatrix<K> &m ) |
---|
| 726 | { |
---|
| 727 | int i,r,c,digits=0,tmp; |
---|
| 728 | |
---|
| 729 | for( i=0; i<m.rows*m.cols; i++ ) |
---|
| 730 | { |
---|
| 731 | tmp = m.a[i].length( ); |
---|
| 732 | |
---|
| 733 | if( tmp > digits ) digits = tmp; |
---|
| 734 | } |
---|
| 735 | |
---|
| 736 | for( r=0; r<m.rows; r++ ) |
---|
| 737 | { |
---|
| 738 | if( m.rows == 1 ) |
---|
| 739 | { |
---|
| 740 | #ifdef KMATRIX_IOSTREAM |
---|
| 741 | s << "<"; |
---|
| 742 | #else |
---|
| 743 | fprintf( stdout,"<" ); |
---|
| 744 | #endif |
---|
| 745 | } |
---|
| 746 | else if( r == 0 ) |
---|
| 747 | { |
---|
| 748 | #ifdef KMATRIX_IOSTREAM |
---|
| 749 | s << "/"; |
---|
| 750 | #else |
---|
| 751 | fprintf( stdout,"/" ); |
---|
| 752 | #endif |
---|
| 753 | } |
---|
| 754 | else if( r == m.rows - 1 ) |
---|
| 755 | { |
---|
| 756 | #ifdef KMATRIX_IOSTREAM |
---|
| 757 | s << "\\"; |
---|
| 758 | #else |
---|
| 759 | fprintf( stdout,"\\" ); |
---|
| 760 | #endif |
---|
| 761 | } |
---|
| 762 | else |
---|
| 763 | { |
---|
| 764 | #ifdef KMATRIX_IOSTREAM |
---|
| 765 | s << "|"; |
---|
| 766 | #else |
---|
| 767 | fprintf( stdout,"|" ); |
---|
| 768 | #endif |
---|
| 769 | } |
---|
| 770 | |
---|
| 771 | for( c=0; c<m.cols; c++ ) |
---|
| 772 | { |
---|
| 773 | #ifdef KMATRIX_IOSTREAM |
---|
| 774 | s << " "; |
---|
| 775 | #else |
---|
| 776 | fprintf( stdout," " ); |
---|
| 777 | #endif |
---|
| 778 | |
---|
| 779 | print_rational( s,digits,m.a[r*m.cols+c] ); |
---|
| 780 | } |
---|
| 781 | |
---|
| 782 | if( m.rows == 1 ) |
---|
| 783 | { |
---|
| 784 | #ifdef KMATRIX_IOSTREAM |
---|
| 785 | s << " >"; |
---|
| 786 | #else |
---|
| 787 | fprintf( stdout," >" ); |
---|
| 788 | #endif |
---|
| 789 | } |
---|
| 790 | else if( r == 0 ) |
---|
| 791 | { |
---|
| 792 | #ifdef KMATRIX_IOSTREAM |
---|
| 793 | s << " \\" << endl; |
---|
| 794 | #else |
---|
| 795 | fprintf( stdout," \\\n" ); |
---|
| 796 | #endif |
---|
| 797 | } |
---|
| 798 | else if( r == m.rows - 1 ) |
---|
| 799 | { |
---|
| 800 | #ifdef KMATRIX_IOSTREAM |
---|
| 801 | s << " /"; |
---|
| 802 | #else |
---|
| 803 | fprintf( stdout," /" ); |
---|
| 804 | #endif |
---|
| 805 | } |
---|
| 806 | else |
---|
| 807 | { |
---|
| 808 | #ifdef KMATRIX_IOSTREAM |
---|
| 809 | s << " |" << endl; |
---|
| 810 | #else |
---|
| 811 | fprintf( stdout," |\n" ); |
---|
| 812 | #endif |
---|
| 813 | } |
---|
| 814 | } |
---|
| 815 | return s; |
---|
| 816 | } |
---|
| 817 | |
---|
| 818 | #endif |
---|
| 819 | |
---|
| 820 | // ---------------------------------------------------------------------------- |
---|
| 821 | // test if the matrix is quadratic |
---|
| 822 | // return value: TRUE or FALSE |
---|
| 823 | // ---------------------------------------------------------------------------- |
---|
| 824 | |
---|
| 825 | template<class K> |
---|
| 826 | int KMatrix<K>::is_quadratic( void ) const |
---|
| 827 | { |
---|
| 828 | return ( rows == cols ? TRUE : FALSE ); |
---|
| 829 | } |
---|
| 830 | |
---|
| 831 | // ---------------------------------------------------------------------------- |
---|
| 832 | // test if the matrix is symmetric |
---|
| 833 | // return value: TRUE or FALSE |
---|
| 834 | // ---------------------------------------------------------------------------- |
---|
| 835 | |
---|
| 836 | template<class K> |
---|
| 837 | int KMatrix<K>::is_symmetric( void ) const |
---|
| 838 | { |
---|
| 839 | if( is_quadratic( ) ) |
---|
| 840 | { |
---|
| 841 | int r,c; |
---|
| 842 | |
---|
| 843 | for( r=1; r<rows; r++ ) |
---|
| 844 | { |
---|
| 845 | for( c=0; c<r; c++ ) |
---|
| 846 | { |
---|
| 847 | if( a[r*cols+c] != a[c*cols+r] ) |
---|
| 848 | { |
---|
| 849 | return FALSE; |
---|
| 850 | } |
---|
| 851 | } |
---|
| 852 | } |
---|
| 853 | return TRUE; |
---|
| 854 | } |
---|
| 855 | else |
---|
| 856 | { |
---|
| 857 | return FALSE; |
---|
| 858 | } |
---|
| 859 | } |
---|
| 860 | |
---|
| 861 | // ---------------------------------------------------------------------------- |
---|
| 862 | // compute the determinant |
---|
| 863 | // return value: the determinant |
---|
| 864 | // ---------------------------------------------------------------------------- |
---|
| 865 | |
---|
| 866 | template<class K> K KMatrix<K>::determinant( void ) const |
---|
| 867 | { |
---|
| 868 | if( !is_quadratic( ) ) |
---|
| 869 | { |
---|
| 870 | return 0; |
---|
| 871 | } |
---|
| 872 | |
---|
| 873 | KMatrix<K> dummy( *this ); |
---|
| 874 | |
---|
| 875 | int r,c,rank = 0; |
---|
| 876 | K g; |
---|
| 877 | K frank,fr; |
---|
| 878 | K det = 1; |
---|
| 879 | |
---|
| 880 | // make sure that the elements of each row have gcd=1 |
---|
| 881 | // this is useful for pivoting |
---|
| 882 | |
---|
| 883 | for( r=0; r<dummy.rows; r++ ) |
---|
| 884 | { |
---|
| 885 | det *= dummy.set_row_primitive( r ); |
---|
| 886 | } |
---|
| 887 | |
---|
| 888 | // search a pivoting element in each column |
---|
| 889 | // perform Gauss elimination |
---|
| 890 | |
---|
| 891 | for( c=0; c<cols && rank<dummy.rows; c++ ) |
---|
| 892 | { |
---|
| 893 | if( ( r = dummy.column_pivot( rank,c )) >= 0 ) |
---|
| 894 | { |
---|
| 895 | det *= dummy.swap_rows( rank,r ); |
---|
| 896 | |
---|
| 897 | for( r=rank+1; r<dummy.rows; r++ ) |
---|
| 898 | { |
---|
| 899 | if( dummy.a[r*cols+c] != (K)0 ) |
---|
| 900 | { |
---|
| 901 | g = gcd( dummy.a[r*cols+c],dummy.a[rank*cols+c] ); |
---|
| 902 | |
---|
| 903 | frank = -dummy.a[r*cols+c]/g; |
---|
| 904 | fr = dummy.a[rank*cols+c]/g; |
---|
| 905 | |
---|
| 906 | det /= dummy.add_rows( rank,r,frank,fr ); |
---|
| 907 | det *= dummy.set_row_primitive( r ); |
---|
| 908 | } |
---|
| 909 | } |
---|
| 910 | rank++; |
---|
| 911 | } |
---|
| 912 | } |
---|
| 913 | |
---|
| 914 | if( rank != dummy.rows ) |
---|
| 915 | { |
---|
| 916 | return 0; |
---|
| 917 | } |
---|
| 918 | |
---|
| 919 | for( r=0; r<dummy.rows; r++ ) |
---|
| 920 | { |
---|
| 921 | det *= dummy.a[r*cols+r]; |
---|
| 922 | } |
---|
| 923 | |
---|
| 924 | return det; |
---|
| 925 | } |
---|
| 926 | |
---|
| 927 | #endif /* KMATRIX_H */ |
---|
| 928 | |
---|
| 929 | // ---------------------------------------------------------------------------- |
---|
| 930 | // kmatrix.h |
---|
| 931 | // end of file |
---|
| 932 | // ---------------------------------------------------------------------------- |
---|