← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/pthreads/misc/ATL_tzsplit.c

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-06-11 15:58:16 UTC
  • mfrom: (1.1.3 upstream)
  • mto: (2.2.21 experimental)
  • mto: This revision was merged to the branch mainline in revision 26.
  • Revision ID: package-import@ubuntu.com-20130611155816-b72z8f621tuhbzn0
Tags: upstream-3.10.1
Import upstream version 3.10.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/pthreads/misc/ATL_tzsplit.c
1
 
/* ---------------------------------------------------------------------
2
 
 *
3
 
 * -- Automatically Tuned Linear Algebra Software (ATLAS)
4
 
 *    (C) Copyright 2000 All Rights Reserved
5
 
 *
6
 
 * -- ATLAS routine -- Version 3.2 -- December 25, 2000
7
 
 *
8
 
 * Author         : Antoine P. Petitet
9
 
 * Originally developed at the University of Tennessee,
10
 
 * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
11
 
 *
12
 
 * ---------------------------------------------------------------------
13
 
 *
14
 
 * -- Copyright notice and Licensing terms:
15
 
 *
16
 
 *  Redistribution  and  use in  source and binary forms, with or without
17
 
 *  modification, are  permitted provided  that the following  conditions
18
 
 *  are met:
19
 
 *
20
 
 * 1. Redistributions  of  source  code  must retain the above copyright
21
 
 *    notice, this list of conditions and the following disclaimer.
22
 
 * 2. Redistributions in binary form must reproduce  the above copyright
23
 
 *    notice,  this list of conditions, and the  following disclaimer in
24
 
 *    the documentation and/or other materials provided with the distri-
25
 
 *    bution.
26
 
 * 3. The name of the University,  the ATLAS group,  or the names of its
27
 
 *    contributors  may not be used to endorse or promote products deri-
28
 
 *    ved from this software without specific written permission.
29
 
 *
30
 
 * -- Disclaimer:
31
 
 *
32
 
 * THIS  SOFTWARE  IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
33
 
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,  INCLUDING,  BUT NOT
34
 
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
35
 
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
36
 
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,  INDIRECT, INCIDENTAL, SPE-
37
 
 * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
38
 
 * TO,  PROCUREMENT  OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
39
 
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
40
 
 * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT  (IN-
41
 
 * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
42
 
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43
 
 *
44
 
 * ---------------------------------------------------------------------
45
 
 */
46
 
/*
47
 
 * Include files
48
 
 */
49
 
#include "atlas_ptmisc.h"
50
 
 
51
 
DIM_TZSPLIT_T ATL_tzsplit
52
 
(
53
 
   const enum ATLAS_UPLO      UPLO,
54
 
   const unsigned int         NT,
55
 
   const int                  M,
56
 
   const int                  N,
57
 
   const int                  K,
58
 
   const int                  NB,
59
 
   unsigned int               * NT1,
60
 
   unsigned int               * NT2,
61
 
   int                        * MNK1,
62
 
   int                        * MNK2
63
 
)
64
 
{
65
 
/*
66
 
 * Purpose
67
 
 * =======
68
 
 *
69
 
 * ATL_tzsplit splits a trapezoidal matrix into two matrices of the same
70
 
 * kind. The initial input matrix sizes are specified by M, N and K:
71
 
 *
72
 
 *              K       N                           N       K
73
 
 *             ______________                 _____________
74
 
 *            |   |          |               |            |\
75
 
 *          M |   |          |             K |            | \
76
 
 *            | __|__________|               |____________|__\
77
 
 *             \  |          |               |            |  |
78
 
 *          K   \ |          |             M |            |  |
79
 
 *               \|__________|               |____________|__|
80
 
 *
81
 
 *           UPLO = AtlasLower               UPLO = AtlasUpper
82
 
 *
83
 
 * This function returns which dimension should be split if any, as well
84
 
 * as the number of threads to be used for each submtrix.
85
 
 *
86
 
 * Arguments
87
 
 * =========
88
 
 *
89
 
 * UPLO    (input)                       const enum ATLAS_UPLO
90
 
 *         On entry,  UPLO  specifies  whether  the  matrix  is upper or
91
 
 *         lower trapezoidal.
92
 
 *
93
 
 * NT      (input)                       const unsigned int
94
 
 *         On entry, NT specifies the initial total number of threads.
95
 
 *
96
 
 * M       (input)                       const int
97
 
 *         On entry, M  specifies  the  number  of  complete rows of the
98
 
 *         trapezoidal matrix.
99
 
 *
100
 
 * N       (input)                       const int
101
 
 *         On entry, N  specifies  the number of complete columns of the
102
 
 *         trapezoidal matrix.
103
 
 *
104
 
 * K       (input)                       const int
105
 
 *         On entry, K specifies the size of the triangular submatrix.
106
 
 *
107
 
 * NB      (input)                       const int
108
 
 *         On entry,  NB  specifies the blocking factor to be used along
109
 
 *         this dimension N.
110
 
 *
111
 
 * NT1     (output)                      unsigned int *
112
 
 *         On exit,  NT1  specifies the number of threads to be used for
113
 
 *         the first part of the problem size MNK1.
114
 
 *
115
 
 * NT2     (output)                      unsigned int *
116
 
 *         On exit,  NT2  specifies the number of threads to be used for
117
 
 *         the second part of the problem size MNK2.
118
 
 *
119
 
 * MNK1    (output)                      int *
120
 
 *         On exit, MNK1 specifies  the length of the problem size to be
121
 
 *         run on the NT1 threads.
122
 
 *
123
 
 * MNK1    (output)                      int *
124
 
 *         On exit, MNK2 specifies  the length of the problem size to be
125
 
 *         run on the NT2 threads.
126
 
 *
127
 
 * ---------------------------------------------------------------------
128
 
 */
129
 
/*
130
 
 * .. Local Variables ..
131
 
 */
132
 
   int                        mnks[8], nts[8];
133
 
   double                     gacol=0.0, garow=0.0, k, k2, m, mnk1, mnk2,
134
 
                              n, nt, nt1, nt2, p, r, rcol=0.0, rrow=0.0,
135
 
                              s, tmnk1, tmnk2, work;
136
 
   DIM_TZSPLIT_T              split = AtlasTzNoSplit;
137
 
   int                        kblks, mblks, nblks, nbm1, ntm1;
138
 
/* ..
139
 
 * .. Executable Statements ..
140
 
 *
141
 
 */
142
 
   nbm1  = NB - 1;
143
 
   mblks = (M+nbm1) / NB; nblks = (N+nbm1) / NB; kblks = (K+nbm1) / NB;
144
 
   r     = (double)(mblks + kblks) * (double)(kblks + nblks);
145
 
 
146
 
   if( ( r < 4.0 ) || ( NT < 2 ) ) return( split );
147
 
 
148
 
   m     = (double)(M);   n     = (double)(N);   k     = (double)(K);
149
 
   k2    = k * k;         nt    = (double)(NT);  ntm1  = NT - 1;
150
 
   work  = 2.0 * ( ( k + m ) * n + ( k * m ) ) + k2;
151
 
 
152
 
   if( ( M > 0 ) || ( K > 0 ) )                       /* Column split */
153
 
   {
154
 
      s     = k + m;
155
 
      p     = k * m;
156
 
      gacol = ( 2.0 * nt * n * s  - work ) / ( nt * s );
157
 
 
158
 
      if( ( N > 1 ) && ( gacol >= 0.0 ) )
159
 
      {                                                    /* split N */
160
 
         mnk1   = ( ( N - (int)(gacol / 2.0) + nbm1 ) / NB ) * NB;
161
 
         mnk1   = Mmin( mnk1, N );
162
 
         mnk1   = Mmax( mnk1, 1 );
163
 
         mnk2   = N - mnk1;
164
 
 
165
 
         tmnk1  = (double)(mnk1);
166
 
         tmnk2  = (double)(mnk2);
167
 
 
168
 
         nt1    = (int)floor( ( ( 2 * s * tmnk1 ) / work ) * nt + 0.5 );
169
 
         nt2    = NT - ( nt1 = Mmin( nt1, ntm1 ) );
170
 
 
171
 
         rcol   = ( ( 2.0 *   s * tmnk1            ) / (double)(nt1) ) -
172
 
                  ( ( 2.0 * ( s * tmnk2 + p ) + k2 ) / (double)(nt2) );
173
 
 
174
 
         if( UPLO == AtlasLower )
175
 
         {
176
 
            mnks[0] = mnk1; mnks[1] = mnk2;
177
 
            nts [0] = nt1;  nts [1] = nt2;
178
 
         }
179
 
         else
180
 
         {
181
 
            mnks[0] = mnk2; mnks[1] = mnk1;
182
 
            nts [0] = nt2;  nts [1] = nt1;
183
 
         }
184
 
      }
185
 
      else
186
 
      {                                                    /* split K */
187
 
         mnk1   = ( ( (int)(s * ( 1.0 - sqrt( 1.0 + ( gacol / s ) ) )) +
188
 
                    nbm1 ) / NB ) * NB;
189
 
         mnk1   = Mmin( mnk1, K );
190
 
         mnk1   = Mmax( mnk1, 1 );
191
 
         mnk2   = K - mnk1;
192
 
 
193
 
         tmnk1  = (double)(mnk1);
194
 
         tmnk2  = (double)(mnk2);
195
 
 
196
 
         r      = work - ( tmnk2 * ( tmnk2 + ( 2.0 * m ) ) );
197
 
         nt1    = (int)floor( ( r / work ) * nt + 0.5 );
198
 
         nt2    = NT - ( nt1 = Mmin( nt1, ntm1 ) );
199
 
 
200
 
         rcol   = 2.0 * ( s * n + tmnk1 * ( m + tmnk2 ) ) +
201
 
                  ( tmnk1 * tmnk1 );
202
 
         rcol  /= (double)(nt1);
203
 
         rcol  -= ( tmnk2 * ( 2.0 * m + tmnk2 ) ) / (double)(nt2);
204
 
 
205
 
         if( UPLO == AtlasLower )
206
 
         {
207
 
            mnks[2] = mnk1; mnks[3] = mnk2;
208
 
            nts [2] = nt1;  nts [3] = nt2;
209
 
         }
210
 
         else
211
 
         {
212
 
            mnks[2] = mnk2; mnks[3] = mnk1;
213
 
            nts [2] = nt2;  nts [3] = nt1;
214
 
         }
215
 
      }
216
 
   }
217
 
                                                         /* Row split */
218
 
   if( ( N > 0 ) || ( K > 0 ) )
219
 
   {
220
 
      s     = k + n;
221
 
      p     = k * n;
222
 
      garow = ( 2.0 * nt * m * s - work ) / ( nt * s );
223
 
 
224
 
      if( ( M > 1 ) && ( garow >= 0.0 ) )
225
 
      {                                                    /* split M */
226
 
         mnk1   = ( ( M - (int)(garow / 2.0) + nbm1 ) / NB ) * NB;
227
 
         mnk1   = Mmin( mnk1, M );
228
 
         mnk1   = Mmax( mnk1, 1 );
229
 
         mnk2   = M - mnk1;
230
 
 
231
 
         tmnk1  = (double)(mnk1);
232
 
         tmnk2  = (double)(mnk2);
233
 
 
234
 
         nt1    = (int)floor( ( ( 2 * s * tmnk1 ) / work ) * nt + 0.5 );
235
 
         nt2    = NT - ( nt1 = Mmin( nt1, ntm1 ) );
236
 
 
237
 
         rrow   = ( ( 2.0 *   s * tmnk1            ) / (double)(nt1) ) -
238
 
                  ( ( 2.0 * ( s * tmnk2 + p ) + k2 ) / (double)(nt2) );
239
 
 
240
 
         if( UPLO == AtlasUpper )
241
 
         {
242
 
            mnks[4] = mnk1; mnks[5] = mnk2;
243
 
            nts [4] = nt1;  nts [5] = nt2;
244
 
         }
245
 
         else
246
 
         {
247
 
            mnks[4] = mnk2; mnks[5] = mnk1;
248
 
            nts [4] = nt2;  nts [5] = nt1;
249
 
         }
250
 
      }
251
 
      else
252
 
      {                                                    /* split K */
253
 
         mnk1   = ( ( (int)(s * ( 1.0 - sqrt( 1.0 + ( garow / s ) ) )) +
254
 
                      nbm1 ) / NB ) * NB;
255
 
         mnk1   = Mmin( mnk1, K );
256
 
         mnk1   = Mmax( mnk1, 1 );
257
 
         mnk2   = K - mnk1;
258
 
 
259
 
         tmnk1  = (double)(mnk1);
260
 
         tmnk2  = (double)(mnk2);
261
 
 
262
 
         r      = work - ( tmnk2 * ( tmnk2 + ( 2.0 * n ) ) );
263
 
         nt1    = (int)floor( ( r / work ) * nt + 0.5 );
264
 
         nt2    = NT - ( nt1 = Mmin( nt1, ntm1 ) );
265
 
 
266
 
         rrow   = 2.0 * ( s * m + tmnk1 * ( n + tmnk2 ) ) +
267
 
                  ( tmnk1 * tmnk1 );
268
 
         rrow  /= (double)(nt1);
269
 
         rrow  -= ( tmnk2 * ( 2.0 * n + tmnk2 ) ) / (double)(nt2);
270
 
 
271
 
         if( UPLO == AtlasUpper )
272
 
         {
273
 
            mnks[6] = mnk1; mnks[7] = mnk2;
274
 
            nts [6] = nt1;  nts [7] = nt2;
275
 
         }
276
 
         else
277
 
         {
278
 
            mnks[6] = mnk2; mnks[7] = mnk1;
279
 
            nts [6] = nt2;  nts [7] = nt1;
280
 
         }
281
 
      }
282
 
   }
283
 
 
284
 
   if( ( ( M > 0 ) || ( K > 0 ) ) && ( ( K > 0 ) || ( N > 0 ) ) )
285
 
   {
286
 
      if( Mabs( rcol ) <= Mabs( rrow ) )              /* Column split */
287
 
      {
288
 
         if( ( N > 1 ) && ( gacol >= 0.0 ) )
289
 
         {
290
 
            split = AtlasTzSplitNcol;
291
 
            *NT1  = nts [0]; *NT2  = nts [1];
292
 
            *MNK1 = mnks[0]; *MNK2 = mnks[1];
293
 
         }
294
 
         else
295
 
         {
296
 
            split = AtlasTzSplitKcol;
297
 
            *NT1  = nts [2]; *NT2  = nts [3];
298
 
            *MNK1 = mnks[2]; *MNK2 = mnks[3];
299
 
         }
300
 
      }
301
 
      else
302
 
      {
303
 
         if( ( M > 1 ) && ( garow >= 0.0 ) )
304
 
         {
305
 
            split = AtlasTzSplitMrow;
306
 
            *NT1  = nts [4]; *NT2  = nts [5];
307
 
            *MNK1 = mnks[4]; *MNK2 = mnks[5];
308
 
         }
309
 
         else
310
 
         {
311
 
            split = AtlasTzSplitKrow;
312
 
            *NT1  = nts [6]; *NT2  = nts [7];
313
 
            *MNK1 = mnks[6]; *MNK2 = mnks[7];
314
 
         }
315
 
      }
316
 
   }
317
 
   else if( ( M > 0 ) || ( K > 0 ) )
318
 
   {
319
 
      if( ( N > 1 ) && ( gacol >= 0.0 ) )
320
 
      {
321
 
         split = AtlasTzSplitNcol;
322
 
         *NT1  = nts [0]; *NT2  = nts [1];
323
 
         *MNK1 = mnks[0]; *MNK2 = mnks[1];
324
 
      }
325
 
      else
326
 
      {
327
 
         split = AtlasTzSplitKcol;
328
 
         *NT1  = nts [2]; *NT2  = nts [3];
329
 
         *MNK1 = mnks[2]; *MNK2 = mnks[3];
330
 
      }
331
 
   }
332
 
   else if( ( N > 0 ) || ( K > 0 ) )
333
 
   {
334
 
      if( ( M > 1 ) && ( garow >= 0.0 ) )
335
 
      {
336
 
         split = AtlasTzSplitMrow;
337
 
         *NT1  = nts [4]; *NT2  = nts [5];
338
 
         *MNK1 = mnks[4]; *MNK2 = mnks[5];
339
 
      }
340
 
      else
341
 
      {
342
 
         split = AtlasTzSplitKrow;
343
 
         *NT1  = nts [6]; *NT2  = nts [7];
344
 
         *MNK1 = mnks[6]; *MNK2 = mnks[7];
345
 
      }
346
 
   }
347
 
   else
348
 
   {
349
 
      split = AtlasTzNoSplit;
350
 
   }
351
 
 
352
 
   return( split );
353
 
/*
354
 
 * End of ATL_tzsplit
355
 
 */
356
 
}