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