source: git/ntl/src/mat_RR.c @ 287cc8

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