source: git/kernel/kmatrix.h @ 576f5b

spielwiese
Last change on this file since 576f5b was 2f6fc61, checked in by Hans Schoenemann <hannes@…>, 14 years ago
compile also without HAVE_FACTORY git-svn-id: file:///usr/local/Singular/svn/trunk@13197 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 23.7 KB
Line 
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
41template<class K> class KMatrix
42{
43
44private:
45
46    K    *a;                                // the entries ot the matrix
47    int rows;                               // number of rows
48    int cols;                               // number of columns
49
50public:
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
107template<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
118template<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
166template<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
177template<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
196template<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
208template<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
233template<class K>
234    inline  KMatrix<K>::KMatrix( )
235{
236    copy_zero( );
237}
238
239// ----------------------------------------------------------------------------
240//  Copy constructor
241// ----------------------------------------------------------------------------
242
243template<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
253template<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
273template<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
289template<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
310template<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
338template<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
353template<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
371template<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
399template<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
422template<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
449template<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
469template<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
489template<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
531template<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
552template<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
598template<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
692template<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
707template<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
725template<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
826template<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
837template<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
867template<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// ----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.