atl_ctrsmk.c

来自「基于Blas CLapck的.用过的人知道是干啥的」· C语言 代码 · 共 1,052 行 · 第 1/2 页

/*
 *             Automatically Tuned Linear Algebra Software v3.8.0
 *                    (C) Copyright 1997 R. Clint Whaley
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *   1. Redistributions of source code must retain the above copyright
 *      notice, this list of conditions and the following disclaimer.
 *   2. Redistributions in binary form must reproduce the above copyright
 *      notice, this list of conditions, and the following disclaimer in the
 *      documentation and/or other materials provided with the distribution.
 *   3. The name of the ATLAS group or the names of its contributers may
 *      not be used to endorse or promote products derived from this
 *      software without specific written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */
#include "atlas_kern3.h"

#include "atlas_prefetch.h"
#ifdef Left_
static void trsmLU_2(const int N, const TYPE *A, TYPE *B, const int ldb)
/*
 * 'Left, 'Upper', with 1 col prefetch, written with all dependencies shown,
 * so that compiler can optimize.
 * A is known to be 2x2, with 1/alpha already applied, diagonals already
 * inverted
 */
{
   const TYPE ar11=*A, ai11=A[1], ar12=A[4], ai12=A[5];
   const TYPE ar22=A[6], ai22=A[7];
   TYPE xr1, xi1, xr2, xi2;
   TYPE t0, p0;
   const int ldb2=ldb+ldb;
   TYPE *bn=B+ldb2;
   const int pfd=ldb2+ldb2;
   int j;

   p0 = B[2];
   for (j=N-1; j; j--) /* stop 1 iteration early to stop prefetch */
   {
      xr2 = p0  ; xi2 = B[3];
      xr1 = *B  ; xi1 = B[1];

      t0 = xr2;
      xr2 = ar22*xr2 - ai22*xi2;
      xi2 = ar22*xi2 + ai22*t0;     p0 = bn[2];

      xr1 -= ar12*xr2 - ai12*xi2;
      xi1 -= ar12*xi2 + ai12*xr2;   ATL_pfl1W(bn+pfd);
      t0 = xr1;
      xr1 = ar11*xr1 - ai11*xi1;
      xi1 = ar11*xi1 + ai11*t0;

      *B   = xr1; B[1] = xi1;
      B[2] = xr2; B[3] = xi2;
      B = bn;
      bn += ldb2;
   }
   xr2 = p0  ; xi2 = B[3];
   xr1 = *B  ; xi1 = B[1];

   t0 = xr2;
   xr2 = ar22*xr2 - ai22*xi2;
   xi2 = ar22*xi2 + ai22*t0;

   xr1 -= ar12*xr2 - ai12*xi2;
   xi1 -= ar12*xi2 + ai12*xr2;
   t0 = xr1;
   xr1 = ar11*xr1 - ai11*xi1;
   xi1 = ar11*xi1 + ai11*t0;

   *B   = xr1; B[1] = xi1;
   B[2] = xr2; B[3] = xi2;
}
static void trsmLU_3(const int N, const TYPE *A, TYPE *B, const int ldb)
/*
 * 'Left, 'Upper', with 1 col prefetch, written with all dependencies shown,
 * so that compiler can optimize.
 * A is known to be 3x3, with 1/alpha already applied, diagonals already
 * inverted
 */
{
   const TYPE ar11=*A, ai11=A[1], ar12=A[6], ai12=A[7], ar13=A[12], ai13=A[13];
   const TYPE ar22=A[ 8], ai22=A[ 9], ar23=A[14], ai23=A[15];
   const TYPE ar33=A[16], ai33=A[17];
   TYPE xr1, xi1, xr2, xi2, xr3, xi3;
   TYPE t0, p0;
   const int ldb2=ldb+ldb;
   TYPE *bn=B+ldb2;
   const int pfd=ldb2+ldb2;
   int j;

   p0 = B[4];
   for (j=N-1; j; j--)
   {
      xr3 = p0  ; xi3 = B[5];
      xr1 = *B  ; xi1 = B[1];
      xr2 = B[2]; xi2 = B[3];

      t0 = xr3;
      xr3 = ar33*xr3 - ai33*xi3;
      xi3 = ar33*xi3 + ai33*t0;

      xr2 -= ar23*xr3 - ai23*xi3;     p0 = bn[4];
      xi2 -= ar23*xi3 + ai23*xr3;
      t0 = xr2;
      xr2 = ar22*xr2 - ai22*xi2;
      xi2 = ar22*xi2 + ai22*t0;

      xr1 -= ar13*xr3 - ai13*xi3;     ATL_pfl1W(bn+pfd);
      xi1 -= ar13*xi3 + ai13*xr3;     ATL_pfl1W(bn+pfd+4);
      xr1 -= ar12*xr2 - ai12*xi2;
      xi1 -= ar12*xi2 + ai12*xr2;
      t0 = xr1;
      xr1 = ar11*xr1 - ai11*xi1;
      xi1 = ar11*xi1 + ai11*t0;

      *B   = xr1; B[1] = xi1;
      B[2] = xr2; B[3] = xi2;
      B[4] = xr3; B[5] = xi3;
      B = bn;
      bn += ldb2;
   }
   xr3 = p0  ; xi3 = B[5];
   xr1 = *B  ; xi1 = B[1];
   xr2 = B[2]; xi2 = B[3];

   t0 = xr3;
   xr3 = ar33*xr3 - ai33*xi3;
   xi3 = ar33*xi3 + ai33*t0;

   xr2 -= ar23*xr3 - ai23*xi3;
   xi2 -= ar23*xi3 + ai23*xr3;
   t0 = xr2;
   xr2 = ar22*xr2 - ai22*xi2;
   xi2 = ar22*xi2 + ai22*t0;

   xr1 -= ar13*xr3 - ai13*xi3;
   xi1 -= ar13*xi3 + ai13*xr3;
   xr1 -= ar12*xr2 - ai12*xi2;
   xi1 -= ar12*xi2 + ai12*xr2;
   t0 = xr1;
   xr1 = ar11*xr1 - ai11*xi1;
   xi1 = ar11*xi1 + ai11*t0;

   *B   = xr1; B[1] = xi1;
   B[2] = xr2; B[3] = xi2;
   B[4] = xr3; B[5] = xi3;
}
static void trsmLU_4(const int N, const TYPE *A, TYPE *B, const int ldb)
/*
 * 'Left, 'Upper', with 1 col prefetch, written with all dependencies shown,
 * so that compiler can optimize.
 * A is known to be 4x4, with 1/alpha already applied, diagonals already
 * inverted
 */
{
   const TYPE ar11=*A, ai11=A[1], ar12=A[8], ai12=A[9], ar13=A[16], ai13=A[17],
              ar14=A[24], ai14=A[25];
   const TYPE ar22=A[10], ai22=A[11], ar23=A[18], ai23=A[19],
              ar24=A[26], ai24=A[27];
   const TYPE ar33=A[20], ai33=A[21], ar34=A[28], ai34=A[29];
   const TYPE ar44=A[30], ai44=A[31];
   TYPE xr1, xi1, xr2, xi2, xr3, xi3, xr4, xi4;
   TYPE t0, p0;
   const int ldb2=ldb+ldb;
   TYPE *bn=B+ldb2;
   const int pfd = ldb2+ldb2;
   int j;

   p0 = B[6];
   for (j=N-1; j; j--)
   {
      xr4 = p0  ; xi4 = B[7];
      xr1 = *B; xi1 = B[1];
      xr3 = B[4]; xi3 = B[5];
      xr2 = B[2]; xi2 = B[3];

      t0 = xr4;
      xr4 = ar44*xr4 - ai44*xi4;
      xi4 = ar44*xi4 + ai44*t0;

      xr3 -= ar34*xr4 - ai34*xi4;
      xi3 -= ar34*xi4 + ai34*xr4;
      t0 = xr3;
      xr3 = ar33*xr3 - ai33*xi3;
      xi3 = ar33*xi3 + ai33*t0;

      xr2 -= ar24*xr4 - ai24*xi4;     p0 = bn[6];
      xi2 -= ar24*xi4 + ai24*xr4;
      xr2 -= ar23*xr3 - ai23*xi3;
      xi2 -= ar23*xi3 + ai23*xr3;
      t0 = xr2;
      xr2 = ar22*xr2 - ai22*xi2;      ATL_pfl1W(bn+pfd);
      xi2 = ar22*xi2 + ai22*t0;       ATL_pfl1W(bn+pfd+4);

      xr1 -= ar14*xr4 - ai14*xi4;
      xi1 -= ar14*xi4 + ai14*xr4;
      xr1 -= ar13*xr3 - ai13*xi3;
      xi1 -= ar13*xi3 + ai13*xr3;
      xr1 -= ar12*xr2 - ai12*xi2;
      xi1 -= ar12*xi2 + ai12*xr2;
      t0 = xr1;
      xr1 = ar11*xr1 - ai11*xi1;
      xi1 = ar11*xi1 + ai11*t0;

      *B   = xr1; B[1] = xi1;
      B[2] = xr2; B[3] = xi2;
      B[4] = xr3; B[5] = xi3;
      B[6] = xr4; B[7] = xi4;
      B = bn;
      bn += ldb2;
   }
   xr4 = p0  ; xi4 = B[7];
   xr1 = *B; xi1 = B[1];
   xr3 = B[4]; xi3 = B[5];
   xr2 = B[2]; xi2 = B[3];

   t0 = xr4;
   xr4 = ar44*xr4 - ai44*xi4;
   xi4 = ar44*xi4 + ai44*t0;

   xr3 -= ar34*xr4 - ai34*xi4;
   xi3 -= ar34*xi4 + ai34*xr4;
   t0 = xr3;
   xr3 = ar33*xr3 - ai33*xi3;
   xi3 = ar33*xi3 + ai33*t0;

   xr2 -= ar24*xr4 - ai24*xi4;
   xi2 -= ar24*xi4 + ai24*xr4;
   xr2 -= ar23*xr3 - ai23*xi3;
   xi2 -= ar23*xi3 + ai23*xr3;
   t0 = xr2;
   xr2 = ar22*xr2 - ai22*xi2;
   xi2 = ar22*xi2 + ai22*t0;

   xr1 -= ar14*xr4 - ai14*xi4;
   xi1 -= ar14*xi4 + ai14*xr4;
   xr1 -= ar13*xr3 - ai13*xi3;
   xi1 -= ar13*xi3 + ai13*xr3;
   xr1 -= ar12*xr2 - ai12*xi2;
   xi1 -= ar12*xi2 + ai12*xr2;
   t0 = xr1;
   xr1 = ar11*xr1 - ai11*xi1;
   xi1 = ar11*xi1 + ai11*t0;

   *B   = xr1; B[1] = xi1;
   B[2] = xr2; B[3] = xi2;
   B[4] = xr3; B[5] = xi3;
   B[6] = xr4; B[7] = xi4;
}
static void trsmLL_2(const int N, const TYPE *A, TYPE *B, const int ldb)
/*
 * 'Left', 'Lower', with 1 column prefetch, written with all dependencies
 * shown, so that the compiler can optimize.
 * A is known to be 2x2, with 1/alpha already applied, diagonals already
 * inverted
 */
{
   const TYPE ar11=*A, ai11=A[1], ar21=A[2], ai21=A[3];
   const TYPE ar22=A[6], ai22=A[7];
   const int ldb2 = ldb+ldb;
   TYPE xr1, xi1, xr2, xi2;
   TYPE t0, p0;
   TYPE *pBn=B+ldb2;
   const int pfd=ldb2+ldb2;
   int j;

   p0 = *B;
   for (j=N-1; j; j--)
   {
      xr1 = p0; xi1 = B[1];
      xr2 = B[2]; xi2 = B[3];

      t0 = xr1;
      xr1 = ar11 * xr1 - ai11 * xi1;
      xi1 = ar11 * xi1 + ai11 * t0;     p0 = *pBn;

      xr2 -= ar21*xr1 - ai21*xi1;
      xi2 -= ar21*xi1 + ai21*xr1;       ATL_pfl1W(pBn+pfd);
      t0 = xr2;
      xr2 = ar22*xr2 - ai22*xi2;
      xi2 = ar22*xi2 + ai22*t0;

      *B   = xr1; B[1] = xi1;
      B[2] = xr2; B[3] = xi2;
      B = pBn;
      pBn += ldb2;
   }
   xr1 = p0; xi1 = B[1];
   xr2 = B[2]; xi2 = B[3];

   t0 = xr1;
   xr1 = ar11 * xr1 - ai11 * xi1;
   xi1 = ar11 * xi1 + ai11 * t0;

   xr2 -= ar21*xr1 - ai21*xi1;
   xi2 -= ar21*xi1 + ai21*xr1;
   t0 = xr2;
   xr2 = ar22*xr2 - ai22*xi2;
   xi2 = ar22*xi2 + ai22*t0;

   *B   = xr1; B[1] = xi1;
   B[2] = xr2; B[3] = xi2;
}
static void trsmLL_3(const int N, const TYPE *A, TYPE *B, const int ldb)
/*
 * 'Left', 'Lower', with 1 column prefetch, written with all dependencies
 * shown, so that the compiler can optimize.
 * A is known to be 3x3, with 1/alpha already applied, diagonals already
 * inverted
 */
{
   const TYPE ar11=*A, ai11=A[1], ar21=A[2], ai21=A[3], ar31=A[4], ai31=A[5];
   const TYPE ar22=A[ 8], ai22=A[ 9], ar32=A[10], ai32=A[11];
   const TYPE ar33=A[16], ai33=A[17];
   const int ldb2 = ldb+ldb;
   TYPE xr1, xi1, xr2, xi2, xr3, xi3;
   TYPE t0, p0;
   TYPE *pBn=B+ldb2;
   const int pfd=ldb2+ldb2;
   int j;

   p0 = *B;
   for (j=N-1; j; j--)
   {
      xr1 = p0; xi1 = B[1];
      xr3 = B[4]; xi3 = B[5];
      xr2 = B[2]; xi2 = B[3];

      t0 = xr1;
      xr1 = ar11 * xr1 - ai11 * xi1;
      xi1 = ar11 * xi1 + ai11 * t0;

      xr2 -= ar21*xr1 - ai21*xi1;
      xi2 -= ar21*xi1 + ai21*xr1;
      t0 = xr2;
      xr2 = ar22*xr2 - ai22*xi2;
      xi2 = ar22*xi2 + ai22*t0;     p0 = *pBn;

      xr3 -= ar31*xr1 - ai31*xi1;
      xi3 -= ar31*xi1 + ai31*xr1;
      xr3 -= ar32*xr2 - ai32*xi2;      ATL_pfl1W(pBn+pfd);
      xi3 -= ar32*xi2 + ai32*xr2;      ATL_pfl1W(pBn+pfd+4);
      t0 = xr3;
      xr3 = ar33*xr3 - ai33*xi3;
      xi3 = ar33*xi3 + ai33*t0;

      *B   = xr1; B[1] = xi1;
      B[2] = xr2; B[3] = xi2;
      B[4] = xr3; B[5] = xi3;
      B = pBn;
      pBn += ldb2;
   }
   xr1 = p0; xi1 = B[1];
   xr3 = B[4]; xi3 = B[5];
   xr2 = B[2]; xi2 = B[3];

   t0 = xr1;
   xr1 = ar11 * xr1 - ai11 * xi1;
   xi1 = ar11 * xi1 + ai11 * t0;

   xr2 -= ar21*xr1 - ai21*xi1;
   xi2 -= ar21*xi1 + ai21*xr1;
   t0 = xr2;
   xr2 = ar22*xr2 - ai22*xi2;
   xi2 = ar22*xi2 + ai22*t0;

   xr3 -= ar31*xr1 - ai31*xi1;
   xi3 -= ar31*xi1 + ai31*xr1;
   xr3 -= ar32*xr2 - ai32*xi2;
   xi3 -= ar32*xi2 + ai32*xr2;
   t0 = xr3;
   xr3 = ar33*xr3 - ai33*xi3;
   xi3 = ar33*xi3 + ai33*t0;

   *B   = xr1; B[1] = xi1;
   B[2] = xr2; B[3] = xi2;
   B[4] = xr3; B[5] = xi3;
}
static void trsmLL_4(const int N, const TYPE *A, TYPE *B, const int ldb)
/*
 * 'Left', 'Lower', with 1 column prefetch, written with all dependencies
 * shown, so that the compiler can optimize.
 * A is known to be 4x4, with 1/alpha already applied, diagonals already
 * inverted
 */
{
   const TYPE ar11=*A, ai11=A[1], ar21=A[2], ai21=A[3], ar31=A[4], ai31=A[5],
              ar41=A[6], ai41=A[7];
   const TYPE ar22=A[10], ai22=A[11], ar32=A[12], ai32=A[13],
              ar42=A[14], ai42=A[15];
   const TYPE ar33=A[20], ai33=A[21], ar43=A[22], ai43=A[23];
   const TYPE ar44=A[30], ai44=A[31];
   const int ldb2 = ldb+ldb;
   TYPE xr1, xi1, xr2, xi2, xr3, xi3, xr4, xi4;
   TYPE t0, p0;
   TYPE *pBn=B+ldb2;
   const int pfd = ldb2+ldb2;
   int j;

   p0 = *B;
   for (j=N-1; j; j--)
   {
      xr1 = p0; xi1 = B[1];
      xr3 = B[4]; xi3 = B[5];
      xr2 = B[2]; xi2 = B[3];
      xr4 = B[6]; xi4 = B[7];

      t0 = xr1;
      xr1 = ar11 * xr1 - ai11 * xi1;
      xi1 = ar11 * xi1 + ai11 * t0;

      xr2 -= ar21*xr1 - ai21*xi1;
      xi2 -= ar21*xi1 + ai21*xr1;
      t0 = xr2;
      xr2 = ar22*xr2 - ai22*xi2;
      xi2 = ar22*xi2 + ai22*t0;

      xr3 -= ar31*xr1 - ai31*xi1;
      xi3 -= ar31*xi1 + ai31*xr1;
      xr3 -= ar32*xr2 - ai32*xi2;     p0 = *pBn;
      xi3 -= ar32*xi2 + ai32*xr2;
      t0 = xr3;
      xr3 = ar33*xr3 - ai33*xi3;
      xi3 = ar33*xi3 + ai33*t0;

      xr4 -= ar41*xr1 - ai41*xi1;     ATL_pfl1W(pBn+pfd);
      xi4 -= ar41*xi1 + ai41*xr1;     ATL_pfl1W(pBn+pfd+4);
      xr4 -= ar42*xr2 - ai42*xi2;
      xi4 -= ar42*xi2 + ai42*xr2;
      xr4 -= ar43*xr3 - ai43*xi3;
      xi4 -= ar43*xi3 + ai43*xr3;
      t0 = xr4;
      xr4 = ar44*xr4 - ai44*xi4;
      xi4 = ar44*xi4 + ai44*t0;
      *B   = xr1; B[1] = xi1;
      B[2] = xr2; B[3] = xi2;
      B[4] = xr3; B[5] = xi3;
      B[6] = xr4; B[7] = xi4;
      B = pBn;
      pBn += ldb2;
   }
   xr1 = p0; xi1 = B[1];
   xr3 = B[4]; xi3 = B[5];
   xr2 = B[2]; xi2 = B[3];
   xr4 = B[6]; xi4 = B[7];

   t0 = xr1;
   xr1 = ar11 * xr1 - ai11 * xi1;
   xi1 = ar11 * xi1 + ai11 * t0;

   xr2 -= ar21*xr1 - ai21*xi1;
   xi2 -= ar21*xi1 + ai21*xr1;
   t0 = xr2;
   xr2 = ar22*xr2 - ai22*xi2;
   xi2 = ar22*xi2 + ai22*t0;

   xr3 -= ar31*xr1 - ai31*xi1;
   xi3 -= ar31*xi1 + ai31*xr1;
   xr3 -= ar32*xr2 - ai32*xi2;
   xi3 -= ar32*xi2 + ai32*xr2;
   t0 = xr3;
   xr3 = ar33*xr3 - ai33*xi3;
   xi3 = ar33*xi3 + ai33*t0;

   xr4 -= ar41*xr1 - ai41*xi1;
   xi4 -= ar41*xi1 + ai41*xr1;
   xr4 -= ar42*xr2 - ai42*xi2;
   xi4 -= ar42*xi2 + ai42*xr2;
   xr4 -= ar43*xr3 - ai43*xi3;
   xi4 -= ar43*xi3 + ai43*xr3;
   t0 = xr4;
   xr4 = ar44*xr4 - ai44*xi4;
   xi4 = ar44*xi4 + ai44*t0;
   *B   = xr1; B[1] = xi1;
   B[2] = xr2; B[3] = xi2;
   B[4] = xr3; B[5] = xi3;
   B[6] = xr4; B[7] = xi4;
}
#endif

#ifdef Right_
static void trsmRU_2(const int M, const TYPE *A, TYPE *B, const int ldb)
/*
 * 'Right', 'Upper', written with all dependencies shown, so that the
 * compiler can optimize.  A is known to be 2x2, with 1/alpha already applied,
 * diagonals already inverted.
 */
{
   const TYPE ar11=*A, ai11=A[1], ar12=A[4], ai12=A[5];
   const TYPE ar22=A[6], ai22=A[7];
   const int ldb2 = ldb+ldb;
   TYPE xr1, xi1, xr2, xi2, t0;
   TYPE *pB0=B, *pB1 = B+ldb2;
   int i;
   #define PFD 8

   for (i=M; i; i--)
   {
      xr1 = *pB0; xr2 = *pB1;
      xi1 = pB0[1]; xi2 = pB1[1];
/*
 *    real sequence:
 *    x1 *= a11;
 *    x2 = (x2 - x1*a12) * a22;
 */
      t0 = xr1;
      xr1 = xr1*ar11 - xi1*ai11;
      xi1 = t0 *ai11 + xi1*ar11;

      xr2 -= xr1*ar12 - xi1*ai12;
      xi2 -= xr1*ai12 + xi1*ar12;     ATL_pfl1W(pB0+PFD);