Changeset 24
- Timestamp:
- 06/19/07 10:08:48 (2 years ago)
- Location:
- include/fflas-ffpack
- Files:
-
- 3 modified
-
ffpack.h (modified) (14 diffs)
-
ffpack_ludivine.inl (modified) (2 diffs)
-
modular-balanced.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
include/fflas-ffpack/ffpack.h
r23 r24 285 285 for (size_t j=0; j<M;++j) 286 286 F.assign(*(X+i*ldx+j), zero); 287 287 288 // X = L^-1 in n^3/3 289 // ftrtri (F, FflasLower, FflasUnit, M, A, lda); 290 // for (size_t i=1; i<M; ++i) 291 // fcopy (F, i, (X+i*ldx), 1, (A+i*lda), 1); 288 292 invL( F, M, A, lda, X, ldx ); 289 293 // X = Q^-1.X is not necessary since Q = Id … … 366 370 * @param lda leading dimension of A 367 371 * @param P the column permutation 368 * @param Q the row permutation372 * @param Qt the transpose of the row permutation Q 369 373 * @param LuTag flag for setting the earling termination if the matrix 370 374 * is singular … … 376 380 const size_t M, const size_t N, 377 381 typename Field::Element * A, const size_t lda, 378 size_t* P, size_t* Q ,382 size_t* P, size_t* Qt, 379 383 const FFPACK_LUDIVINE_TAG LuTag=FfpackLQUP, 380 384 const size_t cutoff=__FFPACK_LUDIVINE_CUTOFF); … … 419 423 const size_t N, typename Field::Element * A, const size_t lda){ 420 424 421 typename Field::Element mone, one; 425 static typename Field::Element one; 426 static typename Field::Element mone; 422 427 F.init(one,1.0); 423 428 F.init(mone,-1.0); 424 if ((N == 1) && (Diag == FflasNonUnit)) 425 F.divin (*A); 426 else{ 429 if (N == 1){ 430 if (Diag == FflasNonUnit) 431 F.invin (*A); 432 } else { 427 433 size_t N1 = N/2; 428 434 size_t N2 = N - N1; … … 467 473 A+N1, lda, A+N1*lda, lda, one, A, lda); 468 474 469 ftrmm (F, FflasRight, FflasLower, N1, N2, one, A + N1*(lda+1), lda, A + N1, lda);470 471 ftrmm (F, FflasLeft, FflasUpper, N2, N1, one, A + N1*(lda+1), lda, A + N1*lda, lda);475 ftrmm (F, FflasRight, FflasLower, FflasNoTrans, FflasNonUnit, N1, N2, one, A + N1*(lda+1), lda, A + N1, lda); 476 477 ftrmm (F, FflasLeft, FflasUpper, FflasNoTrans, FflasUnit, N2, N1, one, A + N1*(lda+1), lda, A + N1*lda, lda); 472 478 473 479 ftrtrm (F, N2, A + N1*(lda+1), lda); … … 479 485 * 480 486 * After the computation A = [ M \ V ] such that AU = C is a column echelon 481 * decomposition of A, with U = P^T [ V ] and C = M + Q [ Ir ] 482 * [ 0 In-r ] [ 0 ] 487 * decomposition of A, with U = P^T [ V (+Ir) ] and C = M //+ Q [ Ir ] 488 * [ 0 In-r ] // [ 0 ] 489 * Qt = Q^T 483 490 */ 484 491 template <class Field> … … 486 493 ColumnEchelonForm (const Field& F, const size_t M, const size_t N, 487 494 typename Field::Element * A, const size_t lda, 488 size_t* P, size_t* Q ){495 size_t* P, size_t* Qt){ 489 496 490 497 typename Field::Element one, mone; … … 493 500 size_t r; 494 501 495 r = LUdivine (F, FflasUnit, M, N, A, lda, P, Q); 502 // write_field(F,cerr<<"A = "<<std::endl,A,M,N,N); 503 504 r = LUdivine (F, FflasUnit, M, N, A, lda, P, Qt); 505 506 //write_field(F,cerr<<"LUP = "<<std::endl,A,M,N,N); 496 507 497 508 ftrtri (F, FflasUpper, FflasUnit, r, A, lda); 498 509 499 ftrmm (F, FflasLeft, FflasUpper, FflasNoTrans, FflasUnit, r, N ,510 ftrmm (F, FflasLeft, FflasUpper, FflasNoTrans, FflasUnit, r, N-r, 500 511 mone, A, lda, A+r, lda); 501 512 … … 509 520 * After the computation A = [ V ] such that AU = R is a reduced column echelon 510 521 * [ M 0 ] 511 * decomposition of A, where U = P^T [ V] and R = Q [ Ir ]522 * decomposition of A, where U = P^T [ V ] and R = Q [ Ir ] 512 523 * [ 0 In-r ] [ M 0 ] 524 * Qt = Q^T 513 525 */ 514 526 … … 517 529 ReducedColumnEchelonForm (const Field& F, const size_t M, const size_t N, 518 530 typename Field::Element * A, const size_t lda, 519 size_t* P, size_t* Q ){531 size_t* P, size_t* Qt){ 520 532 521 533 typename Field::Element one, mone; … … 524 536 size_t r; 525 537 526 r = ColumnEchelonForm (F, M, N, A, lda, P, Q); 527 538 //write_field(F,std::cerr<<"Start : A = "<<std::endl,A,M,N,lda); 539 540 r = ColumnEchelonForm (F, M, N, A, lda, P, Qt); 541 542 //write_field(F,std::cerr<<"After ColEchelon : A = "<<std::endl,A,M,N,lda); 543 544 528 545 // M = Q^T M 529 for (int i= r-1; i>=0; --i){530 if ( Q [i]> (size_t) i ){531 fswap( F, i ,532 A + Q [i]*lda, 1,546 for (int i=0; i<r; ++i){ 547 if ( Qt[i]> (size_t) i ){ 548 fswap( F, i+1, 549 A + Qt[i]*lda, 1, 533 550 A + i*lda, 1 ); 534 551 } 535 552 } 536 553 554 //write_field(F,std::cerr<<"After Permut A = "<<std::endl,A,M,N,lda); 555 537 556 ftrtri (F, FflasLower, FflasNonUnit, r, A, lda); 538 557 539 ftrmm (F, FflasRight, FflasLower, FflasNoTrans, FflasNonUnit, r, N, 540 one, A, lda, A+r, lda); 541 558 //write_field(F,std::cerr<<"After ftrtri : A = "<<std::endl,A,M,N,lda); 559 560 ftrmm (F, FflasRight, FflasLower, FflasNoTrans, FflasNonUnit, M-r, r, 561 one, A, lda, A+r*lda, lda); 562 563 564 //write_field(F,std::cerr<<"After ftrmm : A = "<<std::endl,A,M,N,lda); 565 542 566 ftrtrm (F, r, A, lda); 567 568 //write_field(F,std::cerr<<"After ftrtrm : A = "<<std::endl,A,M,N,lda); 543 569 544 570 return r; … … 643 669 //for (size_t j=0; j<=Q[i]; ++j) 644 670 //F.init( *(L+Q[i]+j*ldl), 0 ); 645 // cerr<<"1 deplacement "<<i<<"<-->"<<Q[i]<<endl;671 //std::cerr<<"1 deplacement "<<i<<"<-->"<<Q[i]<<endl; 646 672 fcopy( F, LM-Q[i]-1, L+Q[i]*(ldl+1)+ldl,ldl, L+(Q[i]+1)*ldl+i, ldl ); 647 673 for ( size_t j=Q[i]*ldl; j<LM*ldl; j+=ldl) … … 650 676 } 651 677 ftrsm( F, Side, FflasLower, FflasNoTrans, FflasUnit, M, N, one, L, ldl , B, ldb); 652 //write_field(F, cerr<<"dans solveLB "<<endl,L,N,N,ldl);678 //write_field(F,std::cerr<<"dans solveLB "<<endl,L,N,N,ldl); 653 679 // Undo the permutation of L 654 680 for (size_t i=0; i<R; ++i){ … … 818 844 // after the multiplication, problem in ftrmm) 819 845 ftrmm (F, FflasRight, FflasLower, FflasNoTrans, FflasUnit, 820 N2, N1, one, X11, ldx, X21, ldx );821 for (size_t i=0; i<N2; ++i)822 for (size_t j=0; j<N1; ++j)823 F.negin(*(X21+i*ldx+j));846 N2, N1, mone, X11, ldx, X21, ldx ); 847 // for (size_t i=0; i<N2; ++i) 848 // for (size_t j=0; j<N1; ++j) 849 // F.negin(*(X21+i*ldx+j)); 824 850 825 851 // X21 = X22^-1 . X21 -
include/fflas-ffpack/ffpack_ludivine.inl
r21 r24 499 499 if (R < Nup){ 500 500 // Permutation of the 0 rows 501 for ( size_t i = Nup, j = R ; i < Nup + R2; ++i, ++j){ 502 fcopy( F, N - j, A + j*(lda + 1), 1, A + i*lda + j, 1); 503 for (typename Field::Element *Ai = A + i*lda + j; 504 Ai != A + i*lda + N; ++Ai) 505 F.assign (*Ai, zero); 506 size_t t = Q[j]; 507 Q[j]=Q[i]; 508 Q[i] = t; 501 if (Diag == FflasNonUnit){ 502 for ( size_t i = Nup, j = R ; i < Nup + R2; ++i, ++j){ 503 fcopy( F, N - j, A + j*(lda + 1), 1, A + i*lda + j, 1); 504 for (typename Field::Element *Ai = A + i*lda + j; 505 Ai != A + i*lda + N; ++Ai) 506 F.assign (*Ai, zero); 507 size_t t = Q[j]; 508 Q[j]=Q[i]; 509 Q[i] = t; 510 } 511 } else { 512 for ( size_t i = Nup, j = R+1 ; i < Nup + R2; ++i, ++j){ 513 fcopy( F, N - j, A + (j-1)*lda + j, 1, A + i*lda + j, 1); 514 for (typename Field::Element *Ai = A + i*lda + j; 515 Ai != A + i*lda + N; ++Ai) 516 F.assign (*Ai, zero); 517 size_t t = Q[j-1]; 518 Q[j-1]=Q[i]; 519 Q[i] = t; 520 } 509 521 } 510 522 } … … 524 536 //--------------------------------------------------------------------- 525 537 526 template <class Field>538 template <class Field> 527 539 size_t 528 540 FFPACK::LUdivine_construct( const Field& F, const FFLAS_DIAG Diag, -
include/fflas-ffpack/modular-balanced.h
r18 r24 323 323 324 324 //tmp = floor (y + 0.5); 325 //tmp = fmod (tmp, modulus);326 tmp = fmod (y, modulus);325 tmp = fmod (tmp, modulus); 326 //tmp = fmod (y, modulus); 327 327 328 328 if ( tmp > half_mod ) return x = tmp - modulus;
