~ubuntu-branches/ubuntu/vivid/atlas/vivid

« back to all changes in this revision

Viewing changes to src/blas/gemm/ATL_mmJITcp.c

  • Committer: Bazaar Package Importer
  • Author(s): Sylvestre Ledru
  • Date: 2009-09-17 23:31:54 UTC
  • mto: (2.2.1 experimental)
  • mto: This revision was merged to the branch mainline in revision 10.
  • Revision ID: james.westby@ubuntu.com-20090917233154-9esw88ub02twbuab
Tags: upstream-3.8.3
ImportĀ upstreamĀ versionĀ 3.8.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/*
 
2
 *             Automatically Tuned Linear Algebra Software v3.8.3
 
3
 *                    (C) Copyright 2007 R. Clint Whaley
 
4
 *
 
5
 * Redistribution and use in source and binary forms, with or without
 
6
 * modification, are permitted provided that the following conditions
 
7
 * are met:
 
8
 *   1. Redistributions of source code must retain the above copyright
 
9
 *      notice, this list of conditions and the following disclaimer.
 
10
 *   2. Redistributions in binary form must reproduce the above copyright
 
11
 *      notice, this list of conditions, and the following disclaimer in the
 
12
 *      documentation and/or other materials provided with the distribution.
 
13
 *   3. The name of the ATLAS group or the names of its contributers may
 
14
 *      not be used to endorse or promote products derived from this
 
15
 *      software without specific written permission.
 
16
 *
 
17
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 
18
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 
19
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 
20
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
 
21
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 
22
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 
23
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 
24
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 
25
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 
26
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 
27
 * POSSIBILITY OF SUCH DAMAGE.
 
28
 *
 
29
 */
 
30
#include "atlas_misc.h"
 
31
#include "atlas_lvl3.h"
 
32
#include <stdlib.h>
 
33
 
 
34
int Mjoin(PATL,mmJITcp)(const enum ATLAS_TRANS TA, const enum ATLAS_TRANS TB,
 
35
                        const int M0, const int N, const int K,
 
36
                        const SCALAR alpha, const TYPE *A, const int lda,
 
37
                        const TYPE *B, const int ldb, const SCALAR beta,
 
38
                        TYPE *C, const int ldc)
 
39
/*
 
40
 * Copy matmul algorithm, copies A and B on-the-fly
 
41
 * If M < 0, allocates only (MB+NB)*KB workspace
 
42
 */
 
43
{
 
44
   void *v=NULL;
 
45
   const TYPE *a=A;
 
46
   TYPE *pA, *pB, *pB0;
 
47
   MAT2BLK2 A2blk, B2blk;
 
48
   NBMM0 NBmm0, NBmm1, pNBmm0;
 
49
   const int M = (M0 >= 0) ? M0 : -M0;
 
50
   int nkblks, nmblks, nnblks, mr, nr, kr, KR, bigK, h, i, j, ZEROC;
 
51
   int incAk, incBk, incAm, incBn, incAW, incAWp, incBW, incBWp, incW;
 
52
 
 
53
/*
 
54
 * If both M and N <= NB, and one of them is not full, call BPP, which
 
55
 * can sometimes avoid doing cleanup forall cases
 
56
 */
 
57
   if (M <= MB && N <= NB && (M != MB || N != NB))
 
58
      return(Mjoin(PATL,mmBPP)(TA, TB, M, N, K, alpha, A, lda, B, ldb,
 
59
                               beta, C, ldc));
 
60
/*
 
61
 * If these workspace increments are 0, we do JIT NBxNB copies instead of
 
62
 * copying entire array/panel.  Don't copy mat if you can't reuse it.
 
63
 */
 
64
   if (M0 > 0)
 
65
   {
 
66
      incAW = (N > NB) ? KB*MB : 0;
 
67
      incBW = (M > NB) ? KB*NB : 0;
 
68
   }
 
69
   else /* allocate in minimal space */
 
70
      incAW = incBW = 0;
 
71
   nmblks = M/MB;
 
72
   nnblks = N/NB;
 
73
   nkblks = K/KB;
 
74
   mr = M - nmblks*MB;
 
75
   nr = N - nnblks*NB;
 
76
   kr = K - nkblks*KB;
 
77
/*
 
78
 * K-loop is special, in that we don't call user cleanup, must explicitly zero,
 
79
 * and K-cleanup is typically slower even for generated kernels.  Therefore,
 
80
 * allow extra leaway for doing extra flops.  Note error is unaffected by
 
81
 * any of these extra flops: K-loop has elts zeroed, and multiplying zeros
 
82
 * and adding in zeros doesn't add to error
 
83
 */
 
84
   KR = (kr && kr+4 >= KB) ? KB : kr;
 
85
   bigK = nkblks*KB+KR;
 
86
   if (incAW)
 
87
   {
 
88
      i = MB*bigK;
 
89
      incAWp = KB*mr;
 
90
   }
 
91
   else
 
92
   {
 
93
      i = MB*KB;
 
94
      incAWp = 0;
 
95
   }
 
96
   if (incBW)
 
97
   {
 
98
      incBWp = KB*nr;
 
99
      incW = bigK*NB;
 
100
      i += N*bigK;
 
101
   }
 
102
   else
 
103
   {
 
104
      incBWp = incW = 0;
 
105
      i += NB*KB;
 
106
   }
 
107
   i *= sizeof(TYPE);
 
108
   if (i <= ATL_MaxMalloc || !(incAW | incBW))
 
109
      v = malloc(ATL_Cachelen+i);
 
110
   if (!v) return(-1);
 
111
   pA = ATL_AlignPtr(v);
 
112
   pB0 = pA + (incAW ? bigK*MB : KB*MB);
 
113
   if (TA == AtlasNoTrans)
 
114
   {
 
115
      A2blk = Mjoin(PATL,gemoveT);
 
116
      incAk = lda*KB;
 
117
      incAm = MB;
 
118
   }
 
119
   else
 
120
   {
 
121
      A2blk = Mjoin(PATL,gemove);
 
122
      incAk = KB;
 
123
      incAm = MB*lda;
 
124
   }
 
125
   if (TB == AtlasNoTrans)
 
126
   {
 
127
      B2blk = Mjoin(PATL,gemove);
 
128
      incBk = KB;
 
129
      incBn = NB*ldb;
 
130
   }
 
131
   else
 
132
   {
 
133
      B2blk = Mjoin(PATL,gemoveT);
 
134
      incBk = ldb*KB;
 
135
      incBn = NB;
 
136
   }
 
137
/*
 
138
 * See what kernel we're calling
 
139
 */
 
140
   if ( SCALAR_IS_ONE(beta) )
 
141
   {
 
142
      NBmm0 = NBmm_b1;
 
143
      pNBmm0 = Mjoin(PATL,pNBmm_b1);
 
144
   }
 
145
   else if ( SCALAR_IS_ZERO(beta) )
 
146
   {
 
147
      NBmm0 = NBmm_b0;
 
148
      pNBmm0 = Mjoin(PATL,pNBmm_b0);
 
149
   }
 
150
   else
 
151
   {
 
152
      NBmm0 = NBmm_bX;
 
153
      pNBmm0 = Mjoin(PATL,pNBmm_bX);
 
154
   }
 
155
   KR = (KR == KB) ? KB : 0;
 
156
   ZEROC = !KR && SCALAR_IS_ZERO(beta);
 
157
 
 
158
   for (i=0; i < nmblks; i++)
 
159
   {
 
160
      a = A+i*incAm;
 
161
      pB = pB0;       /* foreach row-panel of A, start at B's copy space */
 
162
      for (j=nnblks; j; j--)
 
163
      {
 
164
         Mjoin(PATL,mmK)(MB, MB, NB, NB, nkblks, kr, KR, ATL_rone, alpha, beta,
 
165
                         a, lda, incAk, pA, incAW, B, ldb, incBk, pB, incBW,
 
166
                         C, ldc, A2blk, B2blk, NBmm0, NBmm_b1);
 
167
         B += incBn;             /* copy next col panel of B */
 
168
         pB += incW;             /* to next col panel of pB  */
 
169
         a = (incAW ? NULL : a); /* reuse row-panel of A if copied */
 
170
         C += ldc*NB;
 
171
      }
 
172
      if (nr)
 
173
      {
 
174
         if (ZEROC)
 
175
            Mjoin(PATL,gezero)(MB, nr, C, ldc);
 
176
         Mjoin(PATL,mmK)(MB, MB, nr, nr, nkblks, kr, KR, ATL_rone, alpha, beta,
 
177
                         a, lda, incAk, pA, incAW, B, ldb, incBk, pB, incBWp,
 
178
                         C, ldc, A2blk, B2blk, pNBmm0, Mjoin(PATL,pNBmm_b1));
 
179
      }
 
180
      C += MB - nnblks*ldc*NB;
 
181
      if (incBW)
 
182
      {
 
183
         B = NULL;              /* finished copying B */
 
184
         incBn = 0;
 
185
      }
 
186
      else
 
187
         B -= nnblks*incBn;
 
188
   }
 
189
   if (mr)
 
190
   {
 
191
      a = A + nmblks*incAm;
 
192
      pB = pB0;
 
193
      if ( SCALAR_IS_ONE(beta) ) NBmm0 = Mjoin(PATL,pMBmm_b1);
 
194
      else if ( SCALAR_IS_ZERO(beta) ) NBmm0 = Mjoin(PATL,pMBmm_b0);
 
195
      else NBmm0 = Mjoin(PATL,pMBmm_bX);
 
196
      for (j=nnblks; j; j--)
 
197
      {
 
198
         Mjoin(PATL,mmK)(mr, mr, NB, NB, nkblks, kr, KR, ATL_rone, alpha, beta,
 
199
                         a, lda, incAk, pA, incAWp, B, ldb, incBk, pB, incBW,
 
200
                         C, ldc, A2blk, B2blk, NBmm0, Mjoin(PATL,pMBmm_b1));
 
201
         B += incBn;              /* copy next col panel of B */
 
202
         pB += incW;              /* to next col panel of pB  */
 
203
         a = (incAW ? NULL : a);  /* reuse row-panel of A if copied */
 
204
         C += ldc*NB;
 
205
      }
 
206
      if (nr)
 
207
      {
 
208
         if ( SCALAR_IS_ZERO(beta) )
 
209
            Mjoin(PATL,gezero)(mr, nr, C, ldc);
 
210
         Mjoin(PATL,mmK)(mr, mr, nr, nr, nkblks, kr, (incAW | incBW) ? KR:0,
 
211
                         ATL_rone, alpha, beta, a, lda, incAk, pA, incAWp,
 
212
                         B, ldb, incBk, pB, incBWp, C, ldc, A2blk, B2blk,
 
213
                         Mjoin(PATL,pKBmm), Mjoin(PATL,pKBmm));
 
214
      }
 
215
   }
 
216
   free(v);
 
217
   return(0);
 
218
}