Changeset 3 for include/fflas-ffpack/ffpack_charpoly.inl
- Timestamp:
- 03/11/07 15:11:02 (2 years ago)
- Files:
-
- 1 modified
-
include/fflas-ffpack/ffpack_charpoly.inl (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
include/fflas-ffpack/ffpack_charpoly.inl
r1 r3 11 11 template <class Field, class Polynomial> 12 12 std::list<Polynomial>& 13 FFPACK::CharPoly (const Field& F, std::list<Polynomial>& charp, const size_t N,14 15 const enum FFPACK_CHARPOLY_TAG CharpTag){16 switch ( CharpTag) {13 FFPACK::CharPoly (const Field& F, std::list<Polynomial>& charp, const size_t N, 14 typename Field::Element * A, const size_t lda, 15 const enum FFPACK_CHARPOLY_TAG CharpTag){ 16 switch (CharpTag) { 17 17 case FfpackLUK:{ 18 18 typename Field::Element * X = new typename Field::Element[N*(N+1)]; 19 LUKrylov (F, charp, N, A, lda, X, N);19 LUKrylov (F, charp, N, A, lda, X, N); 20 20 delete[] X; 21 21 return charp; 22 22 } 23 23 case FfpackKG:{ 24 return KellerGehrig ( F, charp, N, A, lda);24 return KellerGehrig (F, charp, N, A, lda); 25 25 break; 26 26 } … … 31 31 case FfpackKGFast:{ 32 32 size_t mc, mb, j; 33 if (KGFast ( F, charp, N, A, lda, &mc, &mb, &j)){33 if (KGFast (F, charp, N, A, lda, &mc, &mb, &j)){ 34 34 std::cerr<<"NON GENERIC MATRIX PROVIDED TO KELLER-GEHRIG-FAST"<<std::endl; 35 35 } … … 43 43 case FfpackHybrid:{ 44 44 typename Field::Element * X = new typename Field::Element[N*(N+1)]; 45 LUKrylov_KGFast (F, charp, N, A, lda, X, N);45 LUKrylov_KGFast (F, charp, N, A, lda, X, N); 46 46 delete[] X; 47 return charp; 48 break; 49 } 50 51 case FfpackArithProg:{ 52 size_t attempts=0; 53 bool cont = false; 54 do{ 55 try { 56 CharpolyArithProg (F, charp, N, A, lda, 30); 57 } 58 catch (CharpolyFailed){ 59 if (attempts<2) 60 cont = true; 61 else 62 CharPoly(F, charp, N, A, lda, FfpackLUK); 63 64 } 65 } while (cont); 47 66 return charp; 48 67 break; … … 50 69 default:{ 51 70 typename Field::Element * X = new typename Field::Element[N*(N+1)]; 52 LUKrylov (F, charp, N, A, lda, X, N);71 LUKrylov (F, charp, N, A, lda, X, N); 53 72 delete[] X; 54 73 return charp; … … 60 79 template <class Field, class Polynomial> 61 80 std::list<Polynomial>& 62 FFPACK::LUKrylov (const Field& F, std::list<Polynomial>& charp, const size_t N,63 64 81 FFPACK::LUKrylov (const Field& F, std::list<Polynomial>& charp, const size_t N, 82 typename Field::Element * A, const size_t lda, 83 typename Field::Element * X, const size_t ldx){ 65 84 66 85 typedef typename Field::Element elt; … … 73 92 charp.clear(); 74 93 int nbfac = 0; 75 while ( Ncurr > 0){94 while (Ncurr > 0){ 76 95 size_t P[Ncurr]; 77 Polynomial minP;//=new 78 MinPoly ( F, minP, Ncurr, A, lda, X2, ldx, P);96 Polynomial minP;//=new Polynomial(); 97 MinPoly (F, minP, Ncurr, A, lda, X2, ldx, P); 79 98 int k = minP.size()-1; // degre of minpoly 80 if ( (k==1) && F.isZero( (minP)[0] )){ // minpoly is X99 if ((k==1) && F.isZero ((minP)[0])){ // minpoly is X 81 100 Ai = A; 82 101 int j = Ncurr*Ncurr; 83 while ( j-- && F.isZero(*(Ai++)));84 if ( !j){ // A is 0, CharPoly=X^n102 while (j-- && F.isZero(*(Ai++))); 103 if (!j){ // A is 0, CharPoly=X^n 85 104 minP.resize(Ncurr+1); 86 105 (minP)[1] = zero; … … 90 109 } 91 110 nbfac++; 92 charp.push_front ( minP);93 if ( k==Ncurr){111 charp.push_front (minP); 112 if (k==Ncurr){ 94 113 return charp; 95 114 } … … 97 116 elt * X21 = X2 + k*ldx; 98 117 elt * X22 = X21 + k; 99 // 118 // Compute the n-k last rows of A' = PA^tP^t in X2_ 100 119 // A = A . P^t 101 applyP (F, FflasRight, FflasTrans, Ncurr, 0, k, A, lda, P);120 applyP (F, FflasRight, FflasTrans, Ncurr, 0, k, A, lda, P); 102 121 // Copy X2_ = (A'_2)^t 103 for ( Xi = X21, Ai = A+k; Xi != X21 + Nrest*ldx; Ai++, Xi+=ldx-Ncurr)104 for ( size_t jj=0; jj<Ncurr*lda; jj+=lda)122 for (Xi = X21, Ai = A+k; Xi != X21 + Nrest*ldx; Ai++, Xi+=ldx-Ncurr) 123 for (size_t jj=0; jj<Ncurr*lda; jj+=lda) 105 124 *(Xi++) = *(Ai+jj); 106 125 // A = A . P : Undo the permutation on A 107 applyP (F, FflasRight, FflasNoTrans, Ncurr, 0, k, A, lda, P);108 // X2_ = X2_ . P^t (= ( P A^t P^t )2_)109 applyP (F, FflasRight, FflasTrans, Nrest, 0, k, X21, ldx, P);126 applyP (F, FflasRight, FflasNoTrans, Ncurr, 0, k, A, lda, P); 127 // X2_ = X2_ . P^t (= (P A^t P^t)2_) 128 applyP (F, FflasRight, FflasTrans, Nrest, 0, k, X21, ldx, P); 110 129 // X21 = X21 . S1^-1 111 130 ftrsm(F, FflasRight, FflasUpper, FflasNoTrans, FflasUnit, Nrest, k, 112 one, X2, ldx, X21, ldx); 131 one, X2, ldx, X21, ldx); 113 132 // Creation of the matrix A2 for recurise call 114 for ( Xi = X22,Ai = A;115 116 Xi += (ldx-Nrest), Ai += (lda-Nrest))117 for ( size_t jj=0; jj<Nrest; ++jj)133 for (Xi = X22, Ai = A; 134 Xi != X22 + Nrest*ldx; 135 Xi += (ldx-Nrest), Ai += (lda-Nrest)) 136 for (size_t jj=0; jj<Nrest; ++jj) 118 137 *(Ai++) = *(Xi++); 119 fgemm (F, FflasNoTrans, FflasNoTrans, Nrest, Nrest, k, Mone,120 X21, ldx, X2+k, ldx, one, A, lda 138 fgemm (F, FflasNoTrans, FflasNoTrans, Nrest, Nrest, k, Mone, 139 X21, ldx, X2+k, ldx, one, A, lda); 121 140 X2 = X22; 122 141 Ncurr = Nrest; … … 127 146 template <class Field, class Polynomial> 128 147 std::list<Polynomial>& 129 FFPACK::LUKrylov_KGFast (const Field& F, std::list<Polynomial>& charp, const size_t N,130 131 148 FFPACK::LUKrylov_KGFast (const Field& F, std::list<Polynomial>& charp, const size_t N, 149 typename Field::Element * A, const size_t lda, 150 typename Field::Element * X, const size_t ldx){ 132 151 133 152 typedef typename Field::Element elt; … … 139 158 size_t kg_mc, kg_mb, kg_j; 140 159 141 if (!KGFast ( F, charp, N, A, lda, &kg_mc, &kg_mb, &kg_j))160 if (!KGFast (F, charp, N, A, lda, &kg_mc, &kg_mb, &kg_j)) 142 161 return charp; 143 162 else{// Matrix A is not generic … … 147 166 size_t *P = new size_t[N]; 148 167 149 MinPoly ( F, *minP, N, A, lda, X, ldx, P, FfpackKGF, kg_mc, kg_mb, kg_j);168 MinPoly (F, *minP, N, A, lda, X, ldx, P, FfpackKGF, kg_mc, kg_mb, kg_j); 150 169 size_t k = minP->size()-1; // degre of minpoly 151 if ( (k==1) && F.isZero( (*minP)[0] )){ // minpoly is X170 if ((k==1) && F.isZero ((*minP)[0])){ // minpoly is X 152 171 Ai = A; 153 172 int j = N*N; 154 while ( j-- && F.isZero(*(Ai++)));155 if ( !j){ // A is 0, CharPoly=X^n173 while (j-- && F.isZero(*(Ai++))); 174 if (!j){ // A is 0, CharPoly=X^n 156 175 minP->resize(N+1); 157 176 (*minP)[1] = zero; … … 161 180 } 162 181 163 if ( k==N){182 if (k==N){ 164 183 charp.clear(); 165 184 charp.push_front(*minP); // CharPoly = MinPoly … … 176 195 size_t imax = kg_mc+kg_mb; 177 196 // First Id 178 for ( 197 for (size_t j = 0; j < lambda; ++j){ 179 198 for (size_t i=0; i<imax; ++i) 180 F.assign (*(A+j+i*lda), zero);181 F.assign (*(A+j+imax*lda), one);199 F.assign (*(A+j+i*lda), zero); 200 F.assign (*(A+j+imax*lda), one); 182 201 for (size_t i=imax+1; i<N; ++i) 183 F.assign (*(A+j+i*lda), zero);202 F.assign (*(A+j+i*lda), zero); 184 203 ++imax; 185 204 } 186 205 // Column block B 187 206 for (typename Field::Element* Ai=A; Ai<A+N*lda; Ai+=lda) 188 fcopy ( F, kg_mb, Ai+lambda, 1, Ai+N-kg_mc-kg_mb, 1);207 fcopy (F, kg_mb, Ai+lambda, 1, Ai+N-kg_mc-kg_mb, 1); 189 208 190 209 // Second Id block … … 192 211 for (size_t j = 0; j< kg_j*kg_mc; ++j){ 193 212 for (size_t i = 0; i<imax; ++i) 194 F.assign ( *(A+lambda+kg_mb+j+i*lda), zero);195 F.assign ( *(A+lambda+kg_mb+j+imax*lda), one);213 F.assign (*(A+lambda+kg_mb+j+i*lda), zero); 214 F.assign (*(A+lambda+kg_mb+j+imax*lda), one); 196 215 for (size_t i = imax+1; i<N; ++i) 197 F.assign ( *(A+lambda+kg_mb+j+i*lda), zero);216 F.assign (*(A+lambda+kg_mb+j+i*lda), zero); 198 217 ++imax; 199 218 } 200 219 201 // 220 // Compute the n-k last rows of A' = PA^tP^t in X2_ 202 221 203 222 // A = P . A 204 applyP (F, FflasLeft, FflasNoTrans, N, 0, k,223 applyP (F, FflasLeft, FflasNoTrans, N, 0, k, 205 224 const_cast<typename Field::Element* &>(A), lda, P); 206 225 207 226 // Copy X2_ = (A'2_) 208 for ( Xi = X21, Ai = A+k*lda; Xi != X21 + Nrest*ldx; Ai+=lda-N, Xi+=ldx-N){209 for ( size_t jj=0; jj<N; ++jj){227 for (Xi = X21, Ai = A+k*lda; Xi != X21 + Nrest*ldx; Ai+=lda-N, Xi+=ldx-N){ 228 for (size_t jj=0; jj<N; ++jj){ 210 229 *(Xi++) = *(Ai++); 211 230 } … … 213 232 214 233 // A = P^t . A : Undo the permutation on A 215 applyP (F, FflasLeft, FflasTrans, N, 0, k,234 applyP (F, FflasLeft, FflasTrans, N, 0, k, 216 235 const_cast<typename Field::Element* &>(A), lda, P); 217 236 218 // X2_ = X2_ . P^t (= ( P A P^t )2_)219 applyP (F, FflasRight, FflasTrans, Nrest, 0, k, X21, ldx, P);237 // X2_ = X2_ . P^t (= (P A P^t)2_) 238 applyP (F, FflasRight, FflasTrans, Nrest, 0, k, X21, ldx, P); 220 239 221 240 // X21 = X21 . S1^-1 … … 226 245 elt * A2 = new elt[Nrest*Nrest]; 227 246 228 for ( Xi = X22,A2i = A2;229 230 Xi += (ldx-Nrest)){231 for ( size_t jj=0; jj<Nrest; ++jj){247 for (Xi = X22, A2i = A2; 248 Xi != X22 + Nrest*ldx; 249 Xi += (ldx-Nrest)){ 250 for (size_t jj=0; jj<Nrest; ++jj){ 232 251 *(A2i++) = *(Xi++); 233 252 } 234 253 } 235 fgemm (F, FflasNoTrans, FflasNoTrans, Nrest, Nrest, k, Mone,236 X21, ldx, X+k, ldx, one, A2, Nrest 254 fgemm (F, FflasNoTrans, FflasNoTrans, Nrest, Nrest, k, Mone, 255 X21, ldx, X+k, ldx, one, A2, Nrest); 237 256 238 257 // Recursive call on X22 239 LUKrylov_KGFast ( F, charp, Nrest, A2, Nrest, X22, ldx);240 charp.push_front ( *minP);258 LUKrylov_KGFast (F, charp, Nrest, A2, Nrest, X22, ldx); 259 charp.push_front (*minP); 241 260 delete[] P; 242 261 delete[] A2;
