atl_phk_kmm.c

基于Blas CLapck的.用过的人知道是干啥的

      NBmm_bX(M, N, K, ATL_rone, A, lda, B, ldb, ATL_rnone, C+NBNB, ldc);
      NBmm_b1(M, N, K, ATL_rone, A, lda, B+NBNB, ldb, ATL_rone, C, ldc);
      NBmm_bX(M, N, K, ATL_rone, A+NBNB, lda, B+NBNB, ldb, ATL_rnone,
              C+NBNB, ldc);
      NBmm_b1(M, N, K, ATL_rone, A+NBNB, lda, B, ldb, ATL_rone, C, ldc);
   }
   else if (M != MB)
   {
      if (N == NB && K == KB)
      {
         Mjoin(PATLU,pMBmm_bX)(M, N, K, ATL_rone, A, lda, B, ldb,
                               ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pMBmm_b1)(M, N, K, ATL_rone, A, lda, B+NBNB, ldb,
                               ATL_rone, C, ldc);
         Mjoin(PATLU,pMBmm_bX)(M, N, K, ATL_rone, A+M*K, lda, B+NBNB, ldb,
                               ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pMBmm_b1)(M, N, K, ATL_rone, A+M*K, lda, B, ldb,
                               ATL_rone, C, ldc);
      }
      else
      {
         Mjoin(PATLU,pKBmm)(M, N, K, ATL_rone, A, lda, B, ldb,
                            ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pKBmm)(M, N, K, ATL_rone, A, lda, B+N*K, ldb,
                            ATL_rone, C, ldc);
         Mjoin(PATLU,pKBmm)(M, N, K, ATL_rone, A+M*K, lda, B+N*K, ldb,
                            ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pKBmm)(M, N, K, ATL_rone, A+M*K, lda, B, ldb,
                            ATL_rone, C, ldc);
      }
   }
   else if (N != NB)  /* ib is full */
   {
      if (K == KB)
      {
         Mjoin(PATLU,pNBmm_bX)(M, N, K, ATL_rone, A, lda, B, ldb,
                               ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pNBmm_b1)(M, N, K, ATL_rone, A, lda, B+N*K, ldb,
                               ATL_rone, C, ldc);
         Mjoin(PATLU,pNBmm_bX)(M, N, K, ATL_rone, A+M*K, lda, B+N*K, ldb,
                               ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pNBmm_b1)(M, N, K, ATL_rone, A+M*K, lda, B, ldb,
                               ATL_rone, C, ldc);
      }
      else
      {
         Mjoin(PATLU,pKBmm)(M, N, K, ATL_rone, A, lda, B, ldb,
                            ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pKBmm)(M, N, K, ATL_rone, A, lda, B+N*K, ldb,
                            ATL_rone, C, ldc);
         Mjoin(PATLU,pKBmm)(M, N, K, ATL_rone, A+M*K, lda, B+N*K, ldb,
                            ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pKBmm)(M, N, K, ATL_rone, A+M*K, lda, B, ldb,
                            ATL_rone, C, ldc);
      }
   }
   else  /* ib and jb are full */
   {
         Mjoin(PATLU,pKBmm_bX)(M, N, K, ATL_rone, A, lda, B, ldb,
                               ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pKBmm_b1)(M, N, K, ATL_rone, A, lda, B+N*K, ldb,
                               ATL_rone, C, ldc);
         Mjoin(PATLU,pKBmm_bX)(M, N, K, ATL_rone, A+M*K, lda, B+N*K, ldb,
                               ATL_rnone, C+M*N, ldc);
         Mjoin(PATLU,pKBmm_b1)(M, N, K, ATL_rone, A+M*K, lda, B, ldb,
                               ATL_rone, C, ldc);
   }
}

#define ATL_pmmJIKF      Mjoin(PATL,pmmJIKF)
#define ATL_pmmJIK       Mjoin(PATL,pmmJIK)
#define ATL_prow2blkTF   Mjoin(PATL,prow2blkTF)
#define ATL_pcol2blkF    Mjoin(PATL,pcol2blkF)
#define ATL_pcol2blk     Mjoin(PATL,pcol2blk_a1)
#define ATL_pcol2blk_aX  Mjoin(PATL,pcol2blk_aX)
#define ATL_prow2blkT    Mjoin(PATL,prow2blkT_a1)
#define ATL_prow2blkT_aX Mjoin(PATL,prow2blkT_aX)
#define ATL_pputblk      Mjoin(PATL,pputblk)
#define KBmm            Mjoin(PATL,pKBmm)
#define pNBmm  ATL_gNBmm_b1
#define pNBmm0 ATL_gNBmm_b0
#define ATL_gNBmm ATL_gNBmm_b1
int Mjoin(PATL,phk_kmm)(const enum ATLAS_UPLO UC,
                        const enum PACK_UPLO UA, const enum ATLAS_TRANS TA,
                        const int N, const int K, const SCALAR alpha,
                        const TYPE *A, const int lda, const SCALAR beta,
                        const int CP, TYPE *C, const int ldc)
/*
 * If CP == 0, then C is a general rectangular matrix.
 * If CP != 0, then C is packed matrix.
 * In both cases, UC tells what portion of the matrix to update.
 * Performs a syrk/sprk by calling the gemm kernel directly
 */
{
   const enum PACK_UPLO UC2 = ((CP) ? UC : PackGen);
   const int nKb = ATL_DivByNB(K), kb = K - ATL_MulByNB(nKb);
   const int KK = K - kb;
   const int incK = ATL_MulByNB(K)SHIFT;
   const int ldainc = (UA == AtlasUpper) ? 1 : ((UA == AtlasLower) ? -1 : 0);
   const int ldcinc = (UC2 == AtlasUpper) ? 1 : ((UC2 == AtlasLower) ? -1 : 0);
   int i, j, k, ib, jb, iend;
   void *vC=NULL;
   TYPE *pC, *pA, *pB, *c;
   const int K2 = K + K;
   const TYPE one[2] = {ATL_rone, ATL_rzero};

/*
 * Grab enough space for NBxNB C workspace and all of A, a panel of B,
 * and copy all of A.  If this is too much, we will reduce K in outer routine
 * until malloc succeeds
 */
   i = ATL_Cachelen + ATL_MulBySize(NBNB+(N+NB)*K);
   if (i <= ATL_pkMaxMalloc || K <= NB) vC = malloc(i);
   if (!vC) return(-1);
   pC = ATL_AlignPtr(vC);
   pB = pC + NBNB2;
   pA = pB + incK;
   if (TA == AtlasNoTrans)
      Mjoin(ATL_prow2blkTF,_blk)(NB, N, K, one, A, lda, ldainc, pA);
   else
      Mjoin(Mjoin(PATL,pcol2blkConjF),_blk)(NB, K, N, one, A,lda,ldainc,pA);
/*
 * Loop over column panels of $C$
 */
   for (j=0; j < N; j += NB)
   {
      jb = N-j;
      jb = Mmin(jb, NB);
/*
 *    Conjugate row-panel of A to create col-panel of B
 */
      pA = pC + NBNB2 + incK + j*(K2);
      for (i=0; i < K; i += NB)
      {
         ib = K - i;
         ib = Mmin(ib, NB);
         iend = ib * jb;
         Mjoin(PATLU,cpsc)(iend, -(*alpha), pA, 1, pB, 1);
         Mjoin(PATLU,cpsc)(iend, *alpha, pA+iend, 1, pB+iend, 1);
         pA += iend+iend;
         pB += iend+iend;
      }
      pB = pC + NBNB2;

/*
 *    Loop over row panels of $C$
 */
      iend   = (UC == AtlasUpper) ? j+1 : N;
      for (i=((UC == AtlasUpper) ? 0 : j); i < iend; i += NB)
      {
         pA = pC + NBNB2 + incK + i*K2;
         ib = N-i;
         ib = Mmin(ib, NB);
         if (ib != NB || jb != NB)
            Mjoin(PATL,gezero)(ib, jb, pC, ib);
         if (nKb)
         {
            pNBmm0(ib, jb, NB, ATL_rone, pA, NB, pB, NB, ATL_rzero, pC, ib);
            for (k=NB; k != KK; k += NB)
               pNBmm(ib, jb, NB, ATL_rone, pA+ib*(k SHIFT), NB, pB+jb*(k SHIFT),
                     NB, ATL_rone, pC, ib);
            if (kb)
               ATL_gNBmm(ib, jb, kb, ATL_rone, pA+ib*(KK SHIFT), kb,
                         pB+jb*(KK SHIFT), kb, ATL_rone, pC, ib);
         }
         else if (kb)
         {
            if (ib == NB && jb == NB) Mjoin(PATL,gezero)(ib, jb, pC, ib);
            ATL_gNBmm(ib,jb,kb, ATL_rone, pA, kb, pB, kb, ATL_rzero, pC, ib);
         }
         if (i != j)
            Mjoin(PATL,pputblk_aX)(ib, jb, pC, C+MindexP(UC2,i,j,ldc),
                                   Mpld(UC2,j,ldc), ldcinc, one, beta);
         else
         {
            Mjoin(PATL,pputblk_diag)(ib, jb, pC, UC, C+MindexP(UC2,i,j,ldc),
                                     Mpld(UC2,j,ldc), ldcinc, one, beta);
/*
 *          This is ugly, should specialize putblk_diag to do this.
 *          Zero out imag of main diag; fixes case for non-zero BETA, and
 *          also fixes any computational problem where answer is 0, but
 *          due to roundoff becomes non-zero
 *          NOTE: ib == jb for diag case
 */
            for (k=0; k < jb; k++)
               C[MindexP(UC2,i+k,j+k,ldc)+1] = ATL_rzero;
         }
      }
   }
   free(vC);
   return(0);
}