Changeset 287cc8 in git for ntl/src/mat_lzz_p.c


Ignore:
Timestamp:
Jan 5, 2010, 5:51:13 PM (14 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '2a584933abf2a2d3082034c7586d38bb6de1a30a')
Children:
3c38b3810fd61108b01f123f5a91e13ccff52b20
Parents:
1d43d184dd871d77c1ba8e095d768f22a0fbe92f
Message:
ntl 5.5.2

git-svn-id: file:///usr/local/Singular/svn/trunk@12402 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • ntl/src/mat_lzz_p.c

    r1d43d18 r287cc8  
    1414
    1515
    16  
    17 void add(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B) 
    18 { 
    19    long n = A.NumRows(); 
    20    long m = A.NumCols(); 
    21  
    22    if (B.NumRows() != n || B.NumCols() != m)   
    23       Error("matrix add: dimension mismatch"); 
    24  
    25    X.SetDims(n, m); 
    26  
    27    long i, j; 
    28    for (i = 1; i <= n; i++)   
    29       for (j = 1; j <= m; j++) 
    30          add(X(i,j), A(i,j), B(i,j)); 
    31 } 
    32  
    33 void sub(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B) 
    34 { 
    35    long n = A.NumRows(); 
    36    long m = A.NumCols(); 
    37  
    38    if (B.NumRows() != n || B.NumCols() != m) 
    39       Error("matrix sub: dimension mismatch"); 
    40  
    41    X.SetDims(n, m); 
    42  
    43    long i, j; 
    44    for (i = 1; i <= n; i++) 
    45       for (j = 1; j <= m; j++) 
    46          sub(X(i,j), A(i,j), B(i,j)); 
    47 } 
    48  
     16
     17void add(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B)
     18{
     19   long n = A.NumRows();
     20   long m = A.NumCols();
     21
     22   if (B.NumRows() != n || B.NumCols() != m)
     23      Error("matrix add: dimension mismatch");
     24
     25   X.SetDims(n, m);
     26
     27   long i, j;
     28   for (i = 1; i <= n; i++)
     29      for (j = 1; j <= m; j++)
     30         add(X(i,j), A(i,j), B(i,j));
     31}
     32
     33void sub(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B)
     34{
     35   long n = A.NumRows();
     36   long m = A.NumCols();
     37
     38   if (B.NumRows() != n || B.NumCols() != m)
     39      Error("matrix sub: dimension mismatch");
     40
     41   X.SetDims(n, m);
     42
     43   long i, j;
     44   for (i = 1; i <= n; i++)
     45      for (j = 1; j <= m; j++)
     46         sub(X(i,j), A(i,j), B(i,j));
     47}
     48
    4949
    5050
    … …  
    5555
    5656
    57 static 
    58 void mul_aux(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B) 
    59 { 
    60    long n = A.NumRows(); 
    61    long l = A.NumCols(); 
    62    long m = B.NumCols(); 
    63  
    64    if (l != B.NumRows()) 
    65       Error("matrix mul: dimension mismatch"); 
    66  
    67    X.SetDims(n, m); 
     57static
     58void mul_aux(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B)
     59{
     60   long n = A.NumRows();
     61   long l = A.NumCols();
     62   long m = B.NumCols();
     63
     64   if (l != B.NumRows())
     65      Error("matrix mul: dimension mismatch");
     66
     67   X.SetDims(n, m);
    6868
    6969   if (m > 1) {  // new preconditioning code
    … …  
    8282         for (j = 0; j < m; j++) acc[j] = 0;
    8383
    84          for (k = 0;  k < l; k++) {   
     84         for (k = 0;  k < l; k++) {
    8585            long aa = rep(ap[k]);
    8686            if (aa != 0) {
    … …  
    9292                  T1 = MulModPrecon(rep(bp[j]), aa, p, aapinv);
    9393                  acc[j] = AddMod(acc[j], T1, p);
    94                } 
     94               }
    9595            }
    9696         }
    … …  
    9898         zz_p *xp = X[i].elts();
    9999         for (j = 0; j < m; j++)
    100             xp[j].LoopHole() = acc[j];   
     100            xp[j].LoopHole() = acc[j];
    101101      }
    102102
    … …  
    107107      double pinv = zz_p::ModulusInverse();
    108108
    109       long i, j, k; 
    110       long acc, tmp; 
    111 
    112       for (i = 1; i <= n; i++) { 
    113          for (j = 1; j <= m; j++) { 
    114             acc = 0; 
    115             for(k = 1; k <= l; k++) { 
    116                tmp = MulMod(rep(A(i,k)), rep(B(k,j)), p, pinv); 
    117                acc = AddMod(acc, tmp, p); 
    118             } 
    119             X(i,j).LoopHole() = acc; 
    120          } 
    121       }
    122  
    123    }
    124 } 
    125 
    126 void mul(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B) 
    127 { 
    128    if (&X == &A || &X == &B) { 
    129       mat_zz_p tmp; 
    130       mul_aux(tmp, A, B); 
    131       X = tmp; 
    132    } 
    133    else 
    134       mul_aux(X, A, B); 
    135 } 
     109      long i, j, k;
     110      long acc, tmp;
     111
     112      for (i = 1; i <= n; i++) {
     113         for (j = 1; j <= m; j++) {
     114            acc = 0;
     115            for(k = 1; k <= l; k++) {
     116               tmp = MulMod(rep(A(i,k)), rep(B(k,j)), p, pinv);
     117               acc = AddMod(acc, tmp, p);
     118            }
     119            X(i,j).LoopHole() = acc;
     120         }
     121      }
     122
     123   }
     124}
     125
     126void mul(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B)
     127{
     128   if (&X == &A || &X == &B) {
     129      mat_zz_p tmp;
     130      mul_aux(tmp, A, B);
     131      X = tmp;
     132   }
     133   else
     134      mul_aux(X, A, B);
     135}
    136136
    137137
    … …  
    140140   long l = a.length();
    141141   long m = B.NumCols();
    142  
    143    if (l != B.NumRows()) 
    144       Error("matrix mul: dimension mismatch"); 
    145 
    146    if (m == 0) { 
     142
     143   if (l != B.NumRows())
     144      Error("matrix mul: dimension mismatch");
     145
     146   if (m == 0) {
    147147
    148148      x.SetLength(0);
    149      
     149
    150150   }
    151151   else if (m == 1) {
    … …  
    157157      long k;
    158158
    159       acc = 0; 
    160       for(k = 1; k <= l; k++) { 
    161          tmp = MulMod(rep(a(k)), rep(B(k,1)), p, pinv); 
    162          acc = AddMod(acc, tmp, p); 
    163       } 
     159      acc = 0;
     160      for(k = 1; k <= l; k++) {
     161         tmp = MulMod(rep(a(k)), rep(B(k,1)), p, pinv);
     162         acc = AddMod(acc, tmp, p);
     163      }
    164164
    165165      x.SetLength(1);
    166166      x(1).LoopHole()  = acc;
    167          
     167
    168168   }
    169169   else {  // m > 1.  precondition
    … …  
    194194               acc[j] = AddMod(acc[j], T1, p);
    195195            }
    196          } 
     196         }
    197197      }
    198198
    … …  
    200200      zz_p *xp = x.elts();
    201201      for (j = 0; j < m; j++)
    202          xp[j].LoopHole() = acc[j];   
    203      
    204    }
    205 }
    206 
    207  
     202         xp[j].LoopHole() = acc[j];
     203
     204   }
     205}
     206
     207
    208208void mul_aux(vec_zz_p& x, const mat_zz_p& A, const vec_zz_p& b)
    209209{
    … …  
    258258
    259259         xp[i].LoopHole() = acc;
    260       } 
    261    
    262    }
    263 }
    264  
    265 void mul(vec_zz_p& x, const mat_zz_p& A, const vec_zz_p& b) 
    266 { 
     260      }
     261
     262   }
     263}
     264
     265void mul(vec_zz_p& x, const mat_zz_p& A, const vec_zz_p& b)
     266{
    267267   if (&b == &x || A.position1(x) != -1) {
    268268      vec_zz_p tmp;
    … …  
    273273      mul_aux(x, A, b);
    274274
    275 } 
     275}
    276276
    277277
    … …  
    293293   }
    294294   else {
    295      
     295
    296296      long p = zz_p::modulus();
    297297      double pinv = zz_p::ModulusInverse();
    298298      long bb = rep(b);
    299299      mulmod_precon_t bpinv = PrepMulModPrecon(bb, p, pinv);
    300      
     300
    301301      for (i = 0; i < n; i++) {
    302302         const zz_p *ap = A[i].elts();
    … …  
    315315   b = b_in;
    316316   mul(X, A, b);
    317 } 
    318 
    319 
    320 
    321 
    322      
    323  
    324 void ident(mat_zz_p& X, long n) 
    325 { 
    326    X.SetDims(n, n); 
    327    long i, j; 
    328  
    329    for (i = 1; i <= n; i++) 
    330       for (j = 1; j <= n; j++) 
    331          if (i == j) 
    332             set(X(i, j)); 
    333          else 
    334             clear(X(i, j)); 
    335 } 
     317}
     318
     319
     320
     321
     322
     323
     324void ident(mat_zz_p& X, long n)
     325{
     326   X.SetDims(n, n);
     327   long i, j;
     328
     329   for (i = 1; i <= n; i++)
     330      for (j = 1; j <= n; j++)
     331         if (i == j)
     332            set(X(i, j));
     333         else
     334            clear(X(i, j));
     335}
    336336
    337337
    … …  
    394394
    395395            long T1 = rep(t1);
    396             mulmod_precon_t t1pinv = PrepMulModPrecon(T1, p, pinv); // T1*pinv; 
     396            mulmod_precon_t t1pinv = PrepMulModPrecon(T1, p, pinv); // T1*pinv;
    397397            long T2;
    398398
    … …  
    401401
    402402               T2 = MulModPrecon(rep(*y), T1, p, t1pinv);
    403                x->LoopHole() = AddMod(rep(*x), T2, p); 
     403               x->LoopHole() = AddMod(rep(*x), T2, p);
    404404            }
    405405         }
    … …  
    435435   return 1;
    436436}
    437            
     437
    438438
    439439void transpose(mat_zz_p& X, const mat_zz_p& A)
    … …  
    466466   }
    467467}
    468    
    469 
    470 void solve(zz_p& d, vec_zz_p& X, 
     468
     469
     470void solve(zz_p& d, vec_zz_p& X,
    471471           const mat_zz_p& A, const vec_zz_p& b)
    472472
    … …  
    494494   M.SetDims(n, n+1);
    495495   for (i = 0; i < n; i++) {
    496       for (j = 0; j < n; j++) 
     496      for (j = 0; j < n; j++)
    497497         M[i][j] = A[j][i];
    498498      M[i][n] = b[i];
    … …  
    770770
    771771      D[j] = i;
    772       inv(inverses[j], M[i][j]); 
     772      inv(inverses[j], M[i][j]);
    773773   }
    774774
    … …  
    800800   }
    801801}
    802    
    803 
    804 
    805 
    806 
    807 void diag(mat_zz_p& X, long n, zz_p d) 
    808 { 
    809    X.SetDims(n, n); 
    810    long i, j; 
    811  
    812    for (i = 1; i <= n; i++) 
    813       for (j = 1; j <= n; j++) 
    814          if (i == j) 
    815             X(i, j) = d; 
    816          else 
    817             clear(X(i, j)); 
    818 } 
     802
     803
     804
     805
     806
     807void diag(mat_zz_p& X, long n, zz_p d)
     808{
     809   X.SetDims(n, n);
     810   long i, j;
     811
     812   for (i = 1; i <= n; i++)
     813      for (j = 1; j <= n; j++)
     814         if (i == j)
     815            X(i, j) = d;
     816         else
     817            clear(X(i, j));
     818}
    819819
    820820long IsDiag(const mat_zz_p& A, long n, zz_p d)
Note: See TracChangeset for help on using the changeset viewer.