/ani/mrses

#include <stdio.h>

#include <math.h>

#include <cblas.h>

#include "mrses.h"

int atlas_spotrf2(const int N, MRSESDataType *A, const int lda)

{

   int j;

   MRSESDataType Ajj, *Ac=A, *An=A+lda;

   for (j=0; j != N; j++)

   {

      Ajj = Ac[j] - cblas_sdot(j, Ac, 1, Ac, 1);

      if (Ajj > 0)

      {

         Ac[j] = Ajj = sqrt(Ajj);

         if (j != N-1)

         {

            cblas_sgemv(CblasColMajor, CblasTrans, j, N-j-1, -1,

                       An, lda, Ac, 1, 1, An+j, lda);

            cblas_sscal(N-j-1, 1/Ajj, An+j, lda);

            Ac = An;

            An += lda;

         }

      }

      else

      {

         Ac[j] = Ajj;

         return(j+1);

      }

   }

   return(0);

}

#define TYPE MRSESDataType

#define ATL_rzero 0

#define ATL_rone 1

#define ATL_rnone -1

#define ONE ATL_rone

inline int ATL_potrfU_4(TYPE *A, const short int lda)

{

   TYPE *pA1=A+lda, *pA2=pA1+lda, *pA3=pA2+lda;

   TYPE L11 = *A, L21 = *pA1, L31 = *pA2, L41 = *pA3;

   TYPE L22 = pA1[1], L32 = pA2[1], L42 = pA3[1];

   TYPE L33 = pA2[2], L43 = pA3[2];

   TYPE L44 = pA3[3];

   int iret=0;

   if (L11 > ATL_rzero)

   {

      *A = L11 = sqrt(L11);

      L11 = ATL_rone / L11;

      L21 *= L11;

      L31 *= L11;

      L41 *= L11;

      *pA1 = L21; *pA2 = L31; *pA3 = L41;

      L22 -= L21*L21;

      if (L22 > ATL_rzero)

      {

         pA1[1] = L22 = sqrt(L22);

         L22 = ATL_rone / L22;

         L32 = (L32 - L31*L21) * L22;

         L42 = (L42 - L41*L21) * L22;

         L33 -= L31*L31 + L32*L32;

         pA2[1] = L32; pA3[1] = L42;

         if (L33 > ATL_rzero)

         {

            pA2[2] = L33 = sqrt(L33);

            L43 = (L43 - L41*L31 - L42*L32) / L33;

            L44 -= L41*L41 + L42*L42 + L43*L43;

            pA3[2] = L43;

            if (L44 > ATL_rzero)

            {

               pA3[3] = sqrt(L44);

               return(0);

            }

            else iret=4;

         }

         else iret=3;

      }

      else iret=2;

   }

   else iret=1;

   return(iret);

}

inline int ATL_potrfU_3(TYPE *A, const short int lda)

{

   TYPE *pA1=A+lda, *pA2=pA1+lda;

   register TYPE L11 = *A, L21 = *pA1, L31 = *pA2;

   register TYPE L22=pA1[1], L32=pA2[1];

   register TYPE L33=pA2[2];

   int iret=0;

   if (L11 > ATL_rzero)

   {

      *A = L11 = sqrt(L11);

      L11 = ATL_rone / L11;

      L21 *= L11;

      L31 *= L11;

      *pA1 = L21; *pA2 = L31;

      L22 -= L21*L21;

      if (L22 > ATL_rzero)

      {

         L22 = sqrt(L22);

         L32 = (L32 - L31*L21) / L22;

         L33 -= L31*L31 + L32*L32;

         pA1[1] = L22; pA2[1] = L32;

         if (L33 > ATL_rzero)

         {

            pA2[2] = sqrt(L33);

            return(0);

         }

         else iret=3;

      }

      else iret=2;

   }

   else iret=1;

   return(iret);

}

inline int ATL_potrfU_2(TYPE *A, const short int lda)

{

   TYPE *pA1 = A+lda;

   register TYPE L11=*A, L21=*pA1, L22 = pA1[1];

   if (L11 > ATL_rzero)

   {

      *A = L11 = sqrt(L11);

      *pA1 = L21 = L21 / L11;

      L22 -= L21*L21;

      if (L22 > ATL_rzero)

      {

         pA1[1] = sqrt(L22);

         return(0);

      }

      else return(2);

   }

   else return(1);

}

#define BASE float

#define INDEX short int

inline void

gsl_cblas_strsm_clut_1 (const short int M, const short int N, const float *A, float *B, const short int lda)

{

    register short int i, j, k;

    float Ajj, Bij;

    for (i = 0; i < N; i++) {

      for (j = 0; j < M; j++) {

          Ajj = A[lda * j + j];

          Bij = B[lda * i + j] / Ajj;

	  B[lda * i + j] = Bij;

          for (k = j + 1; k < M; k++) {

            B[lda * i + k] -= A[k * lda + j] * Bij;

          }

      }

    }

}

inline void

gsl_cblas_ssyrk_cut_m11 (const short int N, const short int M, const float *A, float *C, const short int lda) {

  register short int i, j, k;

    for (i = 0; i < N; i++) {

      for (j = 0; j <= i; j++) {

        float temp = 0.0;

        for (k = 0; k < M; k++) {

          temp += A[i * lda + k] * A[j * lda + k];

        }

        C[i * lda + j] -= temp;

      }

    }

}

inline void gsl_spotrf_step(const short int M, const short int N, float *A, const short int lda) {

  register short int i, j, k;

  float *B = A + M*lda;

  float *C = B + M;

  for (i = 0; i < N; i++) {

      for (j = 0; j < M; j++) {

	  float sum = 0;

	  for (k = 0; k < j; k++) {

            sum += A[j * lda + k] *  B[i * lda + k];

          }

	  B[i * lda + j] =  (B[i * lda + j] - sum) / A[j * lda + j];

      }

      for (j = 0; j <= i; j++) {

        float temp = 0.0;

        for (k = 0; k < M; k++) {

          temp += B[i * lda + k] * B[j * lda + k];

	}

        C[i * lda + j] -= temp;

      }

  }

}

int atlas_spotrf_u(const short int N, TYPE *A, const short int lda)

{

   TYPE *Ac, *An;

   short int Nleft, Nright, ierr;

  if (N > 4)

  {

      Nleft = N >> 1;

      Nright = N - Nleft;

      ierr = atlas_spotrf_u(Nleft, A, lda);

      if (!ierr)

      {

        Ac = A + lda * Nleft;

        An = Ac + Nleft; //SHIFT;

        gsl_spotrf_step(Nleft, Nright, A, lda);

/*

//       gsl_cblas_strsm_clut_1(Nleft, Nright, A, Ac, lda);

//	 gsl_cblas_ssyrk_cut_m11(Nright, Nleft, Ac, An, lda);

//         cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans,

//                    CblasNonUnit, Nleft, Nright, ONE, A, lda, Ac, lda);

//         cblas_ssyrk(CblasColMajor, CblasUpper, CblasTrans, Nright, Nleft,

//                  ATL_rnone, Ac, lda, ATL_rone, An, lda);

*/

         ierr = atlas_spotrf_u(Nright, An, lda);

         if (ierr) return(ierr+Nleft);

      }

      else return(ierr);

   }

      else if (N==4) return(ATL_potrfU_4(A, lda));

      else if (N==3) return(ATL_potrfU_3(A, lda));

      else if (N==2) return(ATL_potrfU_2(A, lda));

      else if (N==1)

      {

         if (*A > ATL_rzero) *A = sqrt(*A);

         else return(1);

      }

   return(0);

}