source: git/ntl/src/mat_RR.c @ 2cfffe

spielwiese
Last change on this file since 2cfffe was 2cfffe, checked in by Hans Schönemann <hannes@…>, 21 years ago
This commit was generated by cvs2svn to compensate for changes in r6316, which included commits to RCS files with non-trunk default branches. git-svn-id: file:///usr/local/Singular/svn/trunk@6317 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 11.9 KB
Line 
1
2#include <NTL/mat_RR.h>
3
4#include <NTL/new.h>
5
6NTL_START_IMPL
7
8NTL_matrix_impl(RR,vec_RR,vec_vec_RR,mat_RR)
9NTL_io_matrix_impl(RR,vec_RR,vec_vec_RR,mat_RR)
10NTL_eq_matrix_impl(RR,vec_RR,vec_vec_RR,mat_RR)
11
12
13 
14void add(mat_RR& X, const mat_RR& A, const mat_RR& B) 
15{ 
16   long n = A.NumRows(); 
17   long m = A.NumCols(); 
18 
19   if (B.NumRows() != n || B.NumCols() != m)   
20      Error("matrix add: dimension mismatch"); 
21 
22   X.SetDims(n, m); 
23 
24   long i, j; 
25   for (i = 1; i <= n; i++)   
26      for (j = 1; j <= m; j++) 
27         add(X(i,j), A(i,j), B(i,j)); 
28} 
29 
30void sub(mat_RR& X, const mat_RR& A, const mat_RR& B) 
31{ 
32   long n = A.NumRows(); 
33   long m = A.NumCols(); 
34 
35   if (B.NumRows() != n || B.NumCols() != m) 
36      Error("matrix sub: dimension mismatch"); 
37 
38   X.SetDims(n, m); 
39 
40   long i, j; 
41   for (i = 1; i <= n; i++) 
42      for (j = 1; j <= m; j++) 
43         sub(X(i,j), A(i,j), B(i,j)); 
44} 
45 
46void mul_aux(mat_RR& X, const mat_RR& A, const mat_RR& B) 
47{ 
48   long n = A.NumRows(); 
49   long l = A.NumCols(); 
50   long m = B.NumCols(); 
51 
52   if (l != B.NumRows()) 
53      Error("matrix mul: dimension mismatch"); 
54 
55   X.SetDims(n, m); 
56 
57   long i, j, k; 
58   RR acc, tmp; 
59 
60   for (i = 1; i <= n; i++) { 
61      for (j = 1; j <= m; j++) { 
62         clear(acc); 
63         for(k = 1; k <= l; k++) { 
64            mul(tmp, A(i,k), B(k,j)); 
65            add(acc, acc, tmp); 
66         } 
67         X(i,j) = acc; 
68      } 
69   } 
70} 
71 
72 
73void mul(mat_RR& X, const mat_RR& A, const mat_RR& B) 
74{ 
75   if (&X == &A || &X == &B) { 
76      mat_RR tmp; 
77      mul_aux(tmp, A, B); 
78      X = tmp; 
79   } 
80   else 
81      mul_aux(X, A, B); 
82} 
83 
84 
85static
86void mul_aux(vec_RR& x, const mat_RR& A, const vec_RR& b) 
87{ 
88   long n = A.NumRows(); 
89   long l = A.NumCols(); 
90 
91   if (l != b.length()) 
92      Error("matrix mul: dimension mismatch"); 
93 
94   x.SetLength(n); 
95 
96   long i, k; 
97   RR acc, tmp; 
98 
99   for (i = 1; i <= n; i++) { 
100      clear(acc); 
101      for (k = 1; k <= l; k++) { 
102         mul(tmp, A(i,k), b(k)); 
103         add(acc, acc, tmp); 
104      } 
105      x(i) = acc; 
106   } 
107} 
108 
109 
110void mul(vec_RR& x, const mat_RR& A, const vec_RR& b) 
111{ 
112   if (&b == &x || A.position(b) != -1) {
113      vec_RR tmp;
114      mul_aux(tmp, A, b);
115      x = tmp;
116   }
117   else
118      mul_aux(x, A, b);
119} 
120
121static
122void mul_aux(vec_RR& x, const vec_RR& a, const mat_RR& B) 
123{ 
124   long n = B.NumRows(); 
125   long l = B.NumCols(); 
126 
127   if (n != a.length()) 
128      Error("matrix mul: dimension mismatch"); 
129 
130   x.SetLength(l); 
131 
132   long i, k; 
133   RR acc, tmp; 
134 
135   for (i = 1; i <= l; i++) { 
136      clear(acc); 
137      for (k = 1; k <= n; k++) { 
138         mul(tmp, a(k), B(k,i));
139         add(acc, acc, tmp); 
140      } 
141      x(i) = acc; 
142   } 
143} 
144
145void mul(vec_RR& x, const vec_RR& a, const mat_RR& B)
146{
147   if (&a == &x || B.position(a) != -1) {
148      vec_RR tmp;
149      mul_aux(tmp, a, B);
150      x = tmp;
151   }
152   else
153      mul_aux(x, a, B);
154}
155
156     
157 
158void ident(mat_RR& X, long n) 
159{ 
160   X.SetDims(n, n); 
161   long i, j; 
162 
163   for (i = 1; i <= n; i++) 
164      for (j = 1; j <= n; j++) 
165         if (i == j) 
166            set(X(i, j)); 
167         else 
168            clear(X(i, j)); 
169} 
170
171
172void determinant(RR& d, const mat_RR& M_in)
173{
174   long k, n;
175   long i, j;
176   long pos;
177   RR t1, t2;
178   RR *x, *y;
179
180   n = M_in.NumRows();
181
182   if (M_in.NumCols() != n)
183      Error("determinant: nonsquare matrix");
184
185   if (n == 0) {
186      set(d);
187      return;
188   }
189
190   mat_RR M;
191
192   M = M_in;
193
194
195   RR det;
196   set(det);
197
198   RR maxval;
199
200
201   for (k = 0; k < n; k++) {
202      pos = -1;
203      clear(maxval);
204      for (i = k; i < n; i++) {
205         abs(t1, M[i][k]);
206         if (t1 > maxval) {
207            pos = i;
208            maxval = t1;
209         }
210      }
211
212      if (pos != -1) {
213         if (k != pos) {
214            swap(M[pos], M[k]);
215            negate(det, det);
216         }
217
218         mul(det, det, M[k][k]);
219
220         // make M[k, k] == -1
221
222         inv(t1, M[k][k]);
223         negate(t1, t1);
224         for (j = k+1; j < n; j++) {
225            mul(M[k][j], M[k][j], t1);
226         }
227
228         for (i = k+1; i < n; i++) {
229            // M[i] = M[i] + M[k]*M[i,k]
230
231            t1 = M[i][k];   
232
233            x = M[i].elts() + (k+1);
234            y = M[k].elts() + (k+1);
235
236            for (j = k+1; j < n; j++, x++, y++) {
237               // *x = *x + (*y)*t1
238
239               mul(t2, *y, t1);
240               add(*x, *x, t2);
241            }
242         }
243      }
244      else {
245         clear(d);
246         return;
247      }
248   }
249
250   d = det;
251}
252
253RR determinant(const mat_RR& a)
254   { RR x; determinant(x, a); NTL_OPT_RETURN(RR, x); }
255
256
257long IsIdent(const mat_RR& A, long n)
258{
259   if (A.NumRows() != n || A.NumCols() != n)
260      return 0;
261
262   long i, j;
263
264   for (i = 1; i <= n; i++)
265      for (j = 1; j <= n; j++)
266         if (i != j) {
267            if (!IsZero(A(i, j))) return 0;
268         }
269         else {
270            if (!IsOne(A(i, j))) return 0;
271         }
272
273   return 1;
274}
275           
276
277void transpose(mat_RR& X, const mat_RR& A)
278{
279   long n = A.NumRows();
280   long m = A.NumCols();
281
282   long i, j;
283
284   if (&X == & A) {
285      if (n == m)
286         for (i = 1; i <= n; i++)
287            for (j = i+1; j <= n; j++)
288               swap(X(i, j), X(j, i));
289      else {
290         mat_RR tmp;
291         tmp.SetDims(m, n);
292         for (i = 1; i <= n; i++)
293            for (j = 1; j <= m; j++)
294               tmp(j, i) = A(i, j);
295         X.kill();
296         X = tmp;
297      }
298   }
299   else {
300      X.SetDims(m, n);
301      for (i = 1; i <= n; i++)
302         for (j = 1; j <= m; j++)
303            X(j, i) = A(i, j);
304   }
305}
306   
307
308void solve(RR& d, vec_RR& X, 
309           const mat_RR& A, const vec_RR& b)
310
311{
312   long n = A.NumRows();
313   if (A.NumCols() != n)
314      Error("solve: nonsquare matrix");
315
316   if (b.length() != n)
317      Error("solve: dimension mismatch");
318
319   if (n == 0) {
320      set(d);
321      X.SetLength(0);
322      return;
323   }
324
325   long i, j, k, pos;
326   RR t1, t2;
327   RR *x, *y;
328
329   mat_RR M;
330   M.SetDims(n, n+1);
331
332   for (i = 0; i < n; i++) {
333      for (j = 0; j < n; j++) 
334         M[i][j] = A[j][i];
335      M[i][n] = b[i];
336   }
337
338   RR det;
339   set(det);
340
341   RR maxval;
342
343   for (k = 0; k < n; k++) {
344      pos = -1;
345      clear(maxval);
346      for (i = k; i < n; i++) {
347         abs(t1, M[i][k]);
348         if (t1 > maxval) {
349            pos = i;
350            maxval = t1;
351         }
352      }
353
354      if (pos != -1) {
355         if (k != pos) {
356            swap(M[pos], M[k]);
357            negate(det, det);
358         }
359
360         mul(det, det, M[k][k]);
361
362         // make M[k, k] == -1
363
364         inv(t1, M[k][k]);
365         negate(t1, t1);
366         for (j = k+1; j <= n; j++) {
367            mul(M[k][j], M[k][j], t1);
368         }
369
370         for (i = k+1; i < n; i++) {
371            // M[i] = M[i] + M[k]*M[i,k]
372
373            t1 = M[i][k];   
374
375            x = M[i].elts() + (k+1);
376            y = M[k].elts() + (k+1);
377
378            for (j = k+1; j <= n; j++, x++, y++) {
379               // *x = *x + (*y)*t1
380
381               mul(t2, *y, t1);
382               add(*x, *x, t2);
383            }
384         }
385      }
386      else {
387         clear(d);
388         return;
389      }
390   }
391
392   X.SetLength(n);
393   for (i = n-1; i >= 0; i--) {
394      clear(t1);
395      for (j = i+1; j < n; j++) {
396         mul(t2, X[j], M[i][j]);
397         add(t1, t1, t2);
398      }
399      sub(t1, t1, M[i][n]);
400      X[i] = t1;
401   }
402
403   d = det;
404}
405
406void inv(RR& d, mat_RR& X, const mat_RR& A)
407{
408   long n = A.NumRows();
409   if (A.NumCols() != n)
410      Error("inv: nonsquare matrix");
411
412   if (n == 0) {
413      set(d);
414      X.SetDims(0, 0);
415      return;
416   }
417
418   long i, j, k, pos;
419   RR t1, t2;
420   RR *x, *y;
421
422
423   mat_RR M;
424   M.SetDims(n, 2*n);
425
426   for (i = 0; i < n; i++) {
427      for (j = 0; j < n; j++) {
428         M[i][j] = A[i][j];
429         clear(M[i][n+j]);
430      }
431      set(M[i][n+i]);
432   }
433
434   RR det;
435   set(det);
436
437   RR maxval;
438
439   for (k = 0; k < n; k++) {
440      pos = -1;
441      clear(maxval);
442      for (i = k; i < n; i++) {
443         abs(t1, M[i][k]);
444         if (t1 > maxval) {
445            pos = i;
446            maxval = t1;
447         }
448      }
449
450      if (pos != -1) {
451         if (k != pos) {
452            swap(M[pos], M[k]);
453            negate(det, det);
454         }
455
456         mul(det, det, M[k][k]);
457
458         // make M[k, k] == -1
459
460         inv(t1, M[k][k]);
461         negate(t1, t1);
462         for (j = k+1; j < 2*n; j++) {
463            mul(M[k][j], M[k][j], t1);
464         }
465
466         for (i = k+1; i < n; i++) {
467            // M[i] = M[i] + M[k]*M[i,k]
468
469            t1 = M[i][k];   // this is already reduced
470
471            x = M[i].elts() + (k+1);
472            y = M[k].elts() + (k+1);
473
474            for (j = k+1; j < 2*n; j++, x++, y++) {
475               // *x = *x + (*y)*t1
476
477               mul(t2, *y, t1);
478               add(*x, *x, t2);
479            }
480         }
481      }
482      else {
483         clear(d);
484         return;
485      }
486   }
487
488   X.SetDims(n, n);
489   for (k = 0; k < n; k++) {
490      for (i = n-1; i >= 0; i--) {
491         clear(t1);
492         for (j = i+1; j < n; j++) {
493            mul(t2, X[j][k], M[i][j]);
494            add(t1, t1, t2);
495         }
496         sub(t1, t1, M[i][n+k]);
497         X[i][k] = t1;
498      }
499   }
500
501   d = det;
502}
503
504
505   
506void mul(mat_RR& X, const mat_RR& A, const RR& b_in)
507{
508   RR b = b_in;
509   long n = A.NumRows();
510   long m = A.NumCols();
511
512   X.SetDims(n, m);
513
514   long i, j;
515   for (i = 0; i < n; i++)
516      for (j = 0; j < m; j++)
517         mul(X[i][j], A[i][j], b);
518}
519
520
521void mul(mat_RR& X, const mat_RR& A, double b_in)
522{
523   static RR b;
524   b = b_in;
525   long n = A.NumRows();
526   long m = A.NumCols();
527
528   X.SetDims(n, m);
529
530   long i, j;
531   for (i = 0; i < n; i++)
532      for (j = 0; j < m; j++)
533         mul(X[i][j], A[i][j], b);
534}
535
536void diag(mat_RR& X, long n, const RR& d_in) 
537{ 
538   RR d = d_in;
539   X.SetDims(n, n); 
540   long i, j; 
541 
542   for (i = 1; i <= n; i++) 
543      for (j = 1; j <= n; j++) 
544         if (i == j) 
545            X(i, j) = d; 
546         else 
547            clear(X(i, j)); 
548} 
549
550long IsDiag(const mat_RR& A, long n, const RR& d)
551{
552   if (A.NumRows() != n || A.NumCols() != n)
553      return 0;
554
555   long i, j;
556
557   for (i = 1; i <= n; i++)
558      for (j = 1; j <= n; j++)
559         if (i != j) {
560            if (!IsZero(A(i, j))) return 0;
561         }
562         else {
563            if (A(i, j) != d) return 0;
564         }
565
566   return 1;
567}
568
569void negate(mat_RR& X, const mat_RR& A)
570{
571   long n = A.NumRows();
572   long m = A.NumCols();
573
574
575   X.SetDims(n, m);
576
577   long i, j;
578   for (i = 1; i <= n; i++)
579      for (j = 1; j <= m; j++)
580         negate(X(i,j), A(i,j));
581}
582
583long IsZero(const mat_RR& a)
584{
585   long n = a.NumRows();
586   long i;
587
588   for (i = 0; i < n; i++)
589      if (!IsZero(a[i]))
590         return 0;
591
592   return 1;
593}
594
595void clear(mat_RR& x)
596{
597   long n = x.NumRows();
598   long i;
599   for (i = 0; i < n; i++)
600      clear(x[i]);
601}
602
603
604mat_RR operator+(const mat_RR& a, const mat_RR& b)
605{
606   mat_RR res;
607   add(res, a, b);
608   NTL_OPT_RETURN(mat_RR, res);
609}
610
611mat_RR operator*(const mat_RR& a, const mat_RR& b)
612{
613   mat_RR res;
614   mul_aux(res, a, b);
615   NTL_OPT_RETURN(mat_RR, res);
616}
617
618mat_RR operator-(const mat_RR& a, const mat_RR& b)
619{
620   mat_RR res;
621   sub(res, a, b);
622   NTL_OPT_RETURN(mat_RR, res);
623}
624
625
626mat_RR operator-(const mat_RR& a)
627{
628   mat_RR res;
629   negate(res, a);
630   NTL_OPT_RETURN(mat_RR, res);
631}
632
633
634vec_RR operator*(const mat_RR& a, const vec_RR& b)
635{
636   vec_RR res;
637   mul_aux(res, a, b);
638   NTL_OPT_RETURN(vec_RR, res);
639}
640
641vec_RR operator*(const vec_RR& a, const mat_RR& b)
642{
643   vec_RR res;
644   mul_aux(res, a, b);
645   NTL_OPT_RETURN(vec_RR, res);
646}
647
648
649void inv(mat_RR& X, const mat_RR& A)
650{
651   RR d;
652   inv(d, X, A);
653   if (d == 0) Error("inv: non-invertible matrix");
654}
655
656void power(mat_RR& X, const mat_RR& A, const ZZ& e)
657{
658   if (A.NumRows() != A.NumCols()) Error("power: non-square matrix");
659
660   if (e == 0) {
661      ident(X, A.NumRows());
662      return;
663   }
664
665   mat_RR T1, T2;
666   long i, k;
667
668   k = NumBits(e);
669   T1 = A;
670
671   for (i = k-2; i >= 0; i--) {
672      sqr(T2, T1);
673      if (bit(e, i))
674         mul(T1, T2, A);
675      else
676         T1 = T2;
677   }
678
679   if (e < 0)
680      inv(X, T1);
681   else
682      X = T1;
683}
684
685NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.