Changeset 14

Timestamp:

04/09/07 12:02:27 (2 years ago)

Author:

pernet

Message:

 

Files:

• Makefile (modified) (1 diff)
• include/fflas-ffpack/fflas.h (modified) (1 diff)
• include/fflas-ffpack/ffpack_ludivine.inl (modified) (2 diffs)
• tests/Makefile (modified) (2 diffs)
• tests/test-lqup.C (modified) (5 diffs)

Makefile

@@ -1 +1 @@
 VERSION=1.1.2
 
-FFLAS_FFPACK_DIR=/home/pernet/Logiciels/fflas-ffpack
-LINBOX_DIR=/home/pernet/Logiciels/linbox
+FFLAS_FFPACK_DIR=/home/cpernet/Logiciels/fflas-ffpack
+LINBOX_DIR=/home/cpernet/Logiciels/linbox
 
 

include/fflas-ffpack/fflas.h

@@ -8 +8 @@
  * See COPYING for license information.
  */
+#include <math.h>
 
 #ifndef __FFLAS_H

include/fflas-ffpack/ffpack_ludivine.inl

@@ -37 +37 @@
                         typename Field::Element * A, const size_t lda, size_t*P, 
                         size_t *Q, const enum FFPACK_LUDIVINE_TAG LuTag){
-        static typename Field::Element mone,one;
+        static typename Field::Element mone,one,zero;
         F.init(one,1.0);
+        F.init(zero,0.0);
         F.neg (mone, one);
         size_t MN = MIN(M,N);
@@ -45 +46 @@
         
         for (size_t k = 0; k < MN; ++k){
-                write_field(F,cerr<<"A = "<<endl,A,M,N,lda);
-                cerr<<"Q["<<r<<"] = "<<Q[r]<<endl;
                 size_t p = r;
+                Acurr = A+k*lda+r;
                 while ((p < N) && F.isZero (*(Acurr++)))
                         p++;
-                if (p == N){
-                        Q[r] ++;
-                        Acurr += lda;
-                } else {
+                if (p < N){
                         P[r] = p;
-                        if (r < Q[r])
-                                fcopy (F, N-r, (A + r*(lda+1)), 1, (A+Q[r]*lda+r),1);
+                        if (r < k){
+                                fcopy (F, N-r, (A + r*(lda+1)), 1, (A+k*lda+r),1);
+                                Acurr = A+r+k*lda;
+                                for (size_t i=r; i<N; ++i)
+                                        F.assign(*(Acurr++),zero);
+                        }
+                        
                         fswap (F, M, A+r, lda, A+p, lda);
-                        Acurr += lda +1;
+                        Q[r] = k;
                         r++;
-                        Q[r] = r;
-                }
-                if (Q[r]<M){
-                        write_field(F,cerr<<"A avant  = "<<endl,A,M,N,lda);
-                        ftrsv (F, FflasUpper, FflasTrans, FflasNonUnit, r, A, lda, A+(Q[r])*lda, 1);
-                        write_field(F,cerr<<"A milieu = "<<endl,A,M,N,lda);
-                        fgemv (F, FflasTrans, N-r,r, mone, A+r, lda, A+(Q[r])*lda, 1, one, Acurr, 1);
+                }
+                if (k+1<M){
+                        ftrsv (F, FflasUpper, FflasTrans, FflasNonUnit, r, A, lda, A+(k+1)*lda, 1);
+                        fgemv (F, FflasTrans, r, N-r, mone, A+r, lda, A+(k+1)*lda, 1, one, A+(k+1)*lda+r, 1);
                 } else
                         return r;

tests/Makefile

@@ -4 +4 @@
 # root for the blas library, for ex. /home/foo/ATLAS/lib/Linux_P4SSE2
 #BLASROOT = /home/cpernet/Logiciels/ATLAS_Opteron64/lib/
-BLASROOT = /home/pernet/Logiciels/ATLAS/lib/Linux_PIIISSE2
+BLASROOT = /home/cpernet/Logiciels/GotoBLAS
 
 # ATLAS BLAS users : uncomment these lines:
-CXXFLAGS+=-D__LINBOX_HAVE_CBLAS
-LOADLIBES+=-L${BLASROOT} -lcblas -latlas
+#CXXFLAGS+=-D__LINBOX_HAVE_CBLAS
+#LOADLIBES+=-L${BLASROOT} -lcblas -latlas
 
 # GotoBlas BLAS users : uncomment this line
-#LOADLIBES+=-L${BLASROOT} -lgoto
+LOADLIBES+=-L${BLASROOT} -lgoto
 
 # Other BLAS users, uncomment this line
@@ -17 +17 @@
 
 # architecture parameter for gcc (-march option) for ex. pentium4, athlon ...
-ARCH = -march=pentium3 
+#ARCH = -march=pentium3 
 #ARCH = -march=pentium4
 #ARCH = -march=athlon
 #ARCH = -march=opteron 
+ARCH = -march=opteron -m64 -mtune=k8
 
 #----------------------------------------------------------

tests/test-lqup.C

@@ -20 +20 @@
 #include "Matio.h"
 #include "timer.h"
-#include "fflas-ffpack/modular-balanced.h"
+#include "fflas-ffpack/modular-positive.h"
 #include "fflas-ffpack/ffpack.h"
 
@@ -36 +36 @@
                 exit(-1);
         }
-        Field F(atoi(argv[1]));
+        Field F((unsigned long)atoi(argv[1]));
         Field::Element * A;
         
@@ -63 +63 @@
                 tim.clear();      
                 tim.start();    
-                R = FFPACK::LUdivine (F, diag, m, n, A, n, P, Q,
-                                              FFPACK::FfpackLQUP, cutoff);
+                R = FFPACK::LUdivine_gauss (F, diag, m, n, A, n, P, Q,
+                                            FFPACK::FfpackLQUP);
                 tim.stop();
                 timc+=tim;
         }
         //write_field (F,cerr<<"Result = "<<endl, A, m,n,n);
-        
+
+//      cerr<<"P = [";
+//      for (size_t i=0; i<n; ++i)
+//              cerr<<P[i]<<" ";
+//      cerr<<"]"<<endl;
+//      cerr<<"Q = [";
+//      for (size_t i=0; i<m; ++i)
+//              cerr<<Q[i]<<" ";
+//      cerr<<"]"<<endl;
 #if DEBUG
         Field::Element * L = new Field::Element[m*m];
@@ -84 +92 @@
                         F.assign (*(U + i*n + j), *(A+ i*n+j));
         }
+        for (size_t i=R;i<m; ++i)
+                for (size_t j=0; j<n; ++j)
+                        F.assign(*(U+i*n+j), zero);
         for ( int i=0; i<m; ++i ){
                 int j=0;
-                for (; j< ((i<n)?i:n) ; ++j )
+                for (; j< ((i<R)?i:R) ; ++j )
                         F.assign( *(L + i*m+j), *(A+i*n+j));
                 for (; j<m; ++j )
@@ -92 +103 @@
         }
 
-        FFPACK::applyP( F, FFLAS::FflasRight, FFLAS::FflasNoTrans, m,0,m, L, m, Q);
+//      write_field(F,cerr<<"L = "<<endl,L,m,m,m);
+//      write_field(F,cerr<<"U = "<<endl,U,m,n,n);
+        FFPACK::applyP( F, FFLAS::FflasRight, FFLAS::FflasNoTrans, m,0,R, L, m, Q);
         if (diag == FFLAS::FflasNonUnit)
                 for ( int i=0; i<m; ++i )