Changeset 13

Timestamp:

04/09/07 10:53:44 (2 years ago)

Author:

pernet

Message:

LUdivine Gauss pour experimentations

Location:

include/fflas-ffpack

Files:

• fflas.h (modified) (1 diff)
• ffpack.h (modified) (1 diff)
• ffpack_ludivine.inl (modified) (1 diff)

include/fflas-ffpack/fflas.h

@@ -150 +150 @@
            @brief ftrsv: TRiangular System solve with Vector
            Computes  X <- op(A^-1).X\\
-           size of X is m
+           size of X is N
         */
         template<class Field>

include/fflas-ffpack/ffpack.h

@@ -394 +394 @@
         static size_t 
         LUdivine_small (const Field& F, const enum FFLAS_DIAG Diag,
+                        const size_t M, const size_t N,
+                        typename Field::Element * A, const size_t lda,
+                        size_t* P, size_t* Q, const enum FFPACK_LUDIVINE_TAG LuTag=FfpackLQUP);
+        template <class Field>
+        static size_t 
+        LUdivine_gauss (const Field& F, const enum FFLAS_DIAG Diag,
                         const size_t M, const size_t N,
                         typename Field::Element * A, const size_t lda,

include/fflas-ffpack/ffpack_ludivine.inl

@@ -29 +29 @@
 //              return TURBO (F, M, N, A, lda, P, Q, cutoff);
 // }
+
+
+template<class Field>
+inline size_t
+FFPACK::LUdivine_gauss( const Field& F, const enum FFLAS_DIAG Diag,
+                        const size_t M, const size_t N,         
+                        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;
+        F.init(one,1.0);
+        F.neg (mone, one);
+        size_t MN = MIN(M,N);
+        typename Field::Element * Acurr = A;
+        size_t r = 0;
+        
+        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;
+                while ((p < N) && F.isZero (*(Acurr++)))
+                        p++;
+                if (p == N){
+                        Q[r] ++;
+                        Acurr += lda;
+                } else {
+                        P[r] = p;
+                        if (r < Q[r])
+                                fcopy (F, N-r, (A + r*(lda+1)), 1, (A+Q[r]*lda+r),1);
+                        fswap (F, M, A+r, lda, A+p, lda);
+                        Acurr += lda +1;
+                        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);
+                } else
+                        return r;
+        }
+                        
+        return r;
+}