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 | #include <stdlib.h> |
---|
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 | // ---------------------------------------------------------------------------- |
---|