[Math-atlas-commits] CVS: AtlasBase/Clint atlas-lp.base,1.55,1.56

[R. C. W. <rw...@us...>]

Update of /cvsroot/math-atlas/AtlasBase/Clint
In directory sfp-cvsdas-4.v30.ch3.sourceforge.com:/tmp/cvs-serv22866/Clint

Modified Files:
	atlas-lp.base 
Log Message:


Index: atlas-lp.base
===================================================================
RCS file: /cvsroot/math-atlas/AtlasBase/Clint/atlas-lp.base,v
retrieving revision 1.55
retrieving revision 1.56
diff -C2 -d -r1.55 -r1.56
*** atlas-lp.base	15 Dec 2010 00:56:16 -0000	1.55
--- atlas-lp.base	15 Dec 2010 16:05:46 -0000	1.56
***************
*** 4022,4055 ****
  @beginskip
  int Mjoin(PATL,getf2_r2)
!    (ATL_CINT M, ATL_CINT N, TYPE *A, ATL_CINT lda, int *ipiv)
  {
!    #ifdef TCPLX
!       ATL_CINT lda2 = lda+lda;
!       const TYPE none[2] = {ATL_rnone, ATL_rzero};
!       TYPE alpha[2];
!    #endif
!    ATL_CINT LDA2 = lda2+lda2;
     ATL_CINT MN = Mmin(M,N), N2 = (MN>>1)<<1;
-    int iret=0;
     for (j=0; j < N2; j += 2, Ac += LDA2)
     {
!        if (LU1(M, N, j, A, lda, ipiv))  /* factor column j */
!           iret = j;
  /*
   *     Update column j+1
   */
!        #ifdef TCPLX
!           alpha[0] = -Ac[lda2+j+j];
!           alpha[1] = -Ac[lda2+j+j+1];
!           Mjoin(PATL,axpy)(M-j-1, alpha, Ac+j+j+2, 1, Ac+lda2+j+j+2, 1);
!        #else
!           Mjoin(PATL,axpy)(M-j-1, -Ac[lda+j], Ac+j+1, 1, Ac+lda+j+1, 1);
!        #endif
!        if (LU1(M, N, j+1, A, lda, ipiv))  /* factor column j+1 */
!           iret = j+1;
! /* HERE HERE HERE */
!        my_ger2(M-j-2, N-j-2, none, Ac+((j+2)SHIFT), 1, Ac+((j+2+lda)SHIFT), lda,
!                Ac+LDA2+(2 SHIFT), lda);
     }
  }
  #ifndef TCPLX
--- 4022,4076 ----
  @beginskip
  int Mjoin(PATL,getf2_r2)
!    (ATL_CINT M, ATL_CINT N, TYPE *A0, ATL_CINT lda, int *ipiv)
  {
! #ifdef TCPLX
!    const TYPE none[2] = {ATL_rnone, ATL_rzero};
!    TYPE alpha[2];
!    ATL_CINT lda2 = lda+lda, LDA2 = lda2+lda2;
     ATL_CINT MN = Mmin(M,N), N2 = (MN>>1)<<1;
     for (j=0; j < N2; j += 2, Ac += LDA2)
     {
!       if (LU1(M, N, j, A0, lda, ipiv))  /* factor column j */
!          iret = j;
!    }
  /*
   *     Update column j+1
   */
!        alpha[0] = -Ac[lda2+j+j];
!        alpha[1] = -Ac[lda2+j+j+1];
!        Mjoin(PATL,axpy)(M-j-1, alpha, Ac+j+j+2, 1, Ac+lda2+j+j+2, 1);
! #else
!    TYPE *A = A0;
!    ATL_CINT LDA2 = lda+lda, MN = Mmin(M,N), N2 = (MN>>1)<<1;
!    ATL_CINT incA = lda+lda+2;
!    ATL_INT j;
! 
!    for (j=0; j < N2; j += 2, A += incA)
!    {
!       if (LU1(M, N, j, A0, lda, ipiv))  /* factor column j */
!          iret = j;
! /*
!  *    Update column j+1 and then factor column j+1
!  */
!       Mjoin(PATL,axpy)(M-j-1, -A[lda], A+1, 1, Ac+lda+1, 1);
!       if (LU1(M, N, j+1, A0, lda, ipiv))
!          iret = j+1;
! /*
!  *    Compute 2-wide row-panel of U; L is unit triangular:
!  *    |1   0| |x0| = |b0|  ==> x0 = b0 (no change for A(j,j+1:N))
!  *    |L21 1| |x1|   |b1|  ==> L21*x0 + x1 = b1 ==> x1 = b1 - L21*b0
!  *    So, set alpha = -L21, x = b0, y = x1, and can figure with axpy!
!  */
!       Mjoin(PATL,axpy)(N-j-2, -A[1], A+LDA2, lda, A+LDA2+1, lda);
! /*
!  *    Update trailing matrix with rank-2 update 
!  */
!       my_ger2(M-j-2, N-j-2, none, A+2, 1, A+LDA2, lda,
!               none, A+lda+2, 1, A+LDA2+1, lda, A+incA, lda);
     }
+    if (N-N2)
+       if (LU1(M, N, j, A0, lda, ipiv))  /* factor last column */
+          iret = j;
+    return(iret);
  }
  #ifndef TCPLX