← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/lapack/ATL_getrfC.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/lapack/ATL_getrfC.c
1
1
/*
2
 
 *             Automatically Tuned Linear Algebra Software v3.8.4
 
2
 *             Automatically Tuned Linear Algebra Software v3.10.1
3
3
 *                    (C) Copyright 1999 R. Clint Whaley
4
4
 *
5
5
 * Redistribution and use in source and binary forms, with or without
32
32
#include "atlas_level3.h"
33
33
#include "atlas_level1.h"
34
34
#include "atlas_lapack.h"
 
35
#include "atlas_lamch.h"
 
36
#include Mstr(Mjoin(Mjoin(atlas_,PRE),sysinfo.h))
 
37
#ifdef ATL_USEPTHREADS
 
38
   #include "atlas_pthreads.h"
 
39
   #include "atlas_tcacheedge.h"
 
40
#else
 
41
   #include "atlas_cacheedge.h"
 
42
#endif
 
43
 
35
44
 
36
45
#ifdef TREAL
37
46
   #define ATL_luMmin 2
 
47
   #define ATL_PCAMin 1500
38
48
#else
39
49
   #define ATL_luMmin 1
 
50
   #ifdef SCPLX
 
51
      #define ATL_PCAMin 256
 
52
   #else
 
53
      #define ATL_PCAMin 512
 
54
   #endif
40
55
#endif
41
56
 
42
57
#ifdef TCPLX
45
60
 
46
61
#ifdef TREAL
47
62
 
48
 
#if 0
49
 
static int LU2(const int M, TYPE *A, const int lda, int *ipiv)
50
 
{
51
 
   int ip, iret=0;
52
 
   TYPE *A1 = A + lda;
53
 
   register TYPE t0, t1, scal0, scal1, amax=ATL_rzero;
54
 
   int i, imax=(-1);
55
 
 
56
 
   *ipiv = ip = cblas_iamax(M, A, 1);
57
 
   t0 = A[ip];
58
 
   if (t0 != ATL_rzero)
59
 
   {
60
 
      t1 = A1[ip];
61
 
      A[ip] = *A; A1[ip] = *A1;
62
 
      *A = t0; *A1 = t1;
63
 
      scal0 = ATL_rone / t0; scal1 = -t1;
64
 
      for (i=1; i != M; i++)
65
 
      {
66
 
         t0 = A[i]; t1 = A1[i];
67
 
         t0 *= scal0;
68
 
         t1 += t0 * scal1;
69
 
         A[i] = t0; A1[i] = t1;
70
 
         if (t1 < ATL_rzero) t1 = -t1;
71
 
         if (t1 > amax) { amax = t1; imax = i; }
72
 
      }
73
 
      if (amax != ATL_rzero)
74
 
      {
75
 
         ipiv[1] = imax;
76
 
         t0 = A[imax]; t1 = A1[imax];
77
 
         A[imax] = A[1]; A1[imax] = A1[1];
78
 
         A[1] = t0; A1[1] = t1;
79
 
         cblas_scal(M-2, ATL_rone/t1, A1+2, 1);
80
 
      }
81
 
      else
82
 
      {
83
 
         if (imax != -1) ipiv[1] = imax;
84
 
         else ipiv[1] = 1;
85
 
         iret = 2;
86
 
      }
87
 
   }
88
 
   else iret=1;
89
 
   return(iret);
90
 
}
91
 
#elif 1
92
 
static int LU2(const int M, TYPE *A, const int lda, int *ipiv)
93
 
{
94
 
   int ip, iret=0;
95
 
   TYPE *A1 = A + lda;
96
 
   register TYPE t0, t1, scal0, scal1, amax=ATL_rzero;
97
 
   int i, imax=(-1);
98
 
 
99
 
   *ipiv = ip = cblas_iamax(M, A, 1);
100
 
   t0 = A[ip];
101
 
   if (t0 != ATL_rzero)
102
 
   {
103
 
      t1 = A1[ip];
104
 
      A[ip] = *A; A1[ip] = *A1;
105
 
      *A = t0; *A1 = t1;
106
 
      scal0 = ATL_rone / t0; scal1 = -t1;
107
 
      for (i=1; i != M; i++)
108
 
      {
109
 
         t0 = A[i]; t1 = A1[i];
110
 
         t0 *= scal0;
111
 
         t1 += t0 * scal1;
112
 
         A[i] = t0; A1[i] = t1;
113
 
         if (t1 < ATL_rzero) t1 = -t1;
114
 
         if (t1 > amax) { amax = t1; imax = i; }
115
 
      }
116
 
      if (amax != ATL_rzero)
117
 
      {
118
 
         ipiv[1] = imax;
119
 
         t0 = A[imax]; t1 = A1[imax];
120
 
         A[imax] = A[1]; A1[imax] = A1[1];
121
 
         A[1] = t0; A1[1] = t1;
122
 
         cblas_scal(M-2, ATL_rone/t1, A1+2, 1);
 
63
static int LU2(const int M, TYPE *A, const int lda, int *ipiv)
 
64
/*
 
65
 * factors 2 cols of LU, applies all updated to 2nd call during ininitial
 
66
 * iamax to save time
 
67
 */
 
68
{
 
69
   int ip, iret=0;
 
70
   TYPE *A1 = A + lda;
 
71
   register TYPE t0, t1, scal0, scal1, amax=ATL_rzero;
 
72
   int i, imax=(-1);
 
73
 
 
74
   *ipiv = ip = cblas_iamax(M, A, 1);
 
75
   t0 = A[ip];
 
76
   if (t0 != ATL_rzero)
 
77
   {
 
78
      if (Mabs(t0) >= ATL_laSAFMIN)
 
79
      {
 
80
         t1 = A1[ip];
 
81
         A[ip] = *A; A1[ip] = *A1;
 
82
         *A = t0; *A1 = t1;
 
83
         scal0 = ATL_rone / t0; scal1 = -t1;
 
84
         for (i=1; i != M; i++)
 
85
         {
 
86
            t0 = A[i]; t1 = A1[i];
 
87
            t0 *= scal0;
 
88
            t1 += t0 * scal1;
 
89
            A[i] = t0; A1[i] = t1;
 
90
            if (t1 < ATL_rzero) t1 = -t1;
 
91
            if (t1 > amax) { amax = t1; imax = i; }
 
92
         }
 
93
      }
 
94
      else  /* can't safely invert pivot, so must actually divide */
 
95
      {
 
96
         t1 = A1[ip];
 
97
         A[ip] = *A; A1[ip] = *A1;
 
98
         *A = t0; *A1 = t1;
 
99
         scal0 = t0; scal1 = -t1;
 
100
         for (i=1; i != M; i++)
 
101
         {
 
102
            t0 = A[i]; t1 = A1[i];
 
103
            t0 /= scal0;
 
104
            t1 += t0 * scal1;
 
105
            A[i] = t0; A1[i] = t1;
 
106
            if (t1 < ATL_rzero) t1 = -t1;
 
107
            if (t1 > amax) { amax = t1; imax = i; }
 
108
         }
 
109
      }
 
110
      if (amax != ATL_rzero)
 
111
      {
 
112
         ipiv[1] = imax;
 
113
         t0 = A[imax]; t1 = A1[imax];
 
114
         A[imax] = A[1]; A1[imax] = A1[1];
 
115
         A[1] = t0; A1[1] = t1;
 
116
         if (Mabs(t1) >= ATL_laSAFMIN)
 
117
            cblas_scal(M-2, ATL_rone/t1, A1+2, 1);
 
118
         else
 
119
            for (i=2; i < M; i++)
 
120
               A1[i] /= t1;
123
121
      }
124
122
      else
125
123
      {
139
137
         t0 = A[imax]; t1 = A1[imax];
140
138
         A[imax] = A[1]; A1[imax] = A1[1];
141
139
         A[1] = t0; A1[1] = t1;
142
 
         cblas_scal(M-2, ATL_rone/t1, A1+2, 1);
 
140
         if (Mabs(t1) >= ATL_laSAFMIN)
 
141
            cblas_scal(M-2, ATL_rone/t1, A1+2, 1);
 
142
         else
 
143
            for (i=2; i < M; i++)
 
144
               A1[i] /= t1;
143
145
      }
144
146
      else
145
147
      {
149
151
   }
150
152
   return(iret);
151
153
}
152
 
#else
153
 
static int LU2(const int M, TYPE *A, const int lda, int *ipiv)
154
 
{
155
 
   int ip, iret=0;
156
 
   TYPE *A1 = A + lda;
157
 
   register TYPE t0, t1, scal0, scal1, amax=ATL_rzero;
158
 
   int i, imax=(-1);
159
 
 
160
 
   *ipiv = ip = cblas_iamax(M, A, 1);
161
 
   t0 = A[ip];
162
 
   if (t0 != ATL_rzero)
163
 
   {
164
 
      t1 = A1[ip];
165
 
      A[ip] = *A; A1[ip] = *A1;
166
 
      *A = t0; *A1 = t1;
167
 
      scal0 = ATL_rone / t0; scal1 = -t1;
168
 
   }
169
 
   else
170
 
   {
171
 
      iret = 1;
172
 
      scal0 = ATL_rone; scal1 = -(*A1);
173
 
   }
174
 
   for (i=1; i != M; i++)
175
 
   {
176
 
      t0 = A[i]; t1 = A1[i];
177
 
      t0 *= scal0;
178
 
      t1 += t0 * scal1;
179
 
      A[i] = t0; A1[i] = t1;
180
 
      if (t1 < ATL_rzero) t1 = -t1;
181
 
      if (t1 > amax) { amax = t1; imax = i; }
182
 
   }
183
 
   if (amax != ATL_rzero)
184
 
   {
185
 
      ipiv[1] = imax;
186
 
      t0 = A[imax]; t1 = A1[imax];
187
 
      A[imax] = A[1]; A1[imax] = A1[1];
188
 
      A[1] = t0; A1[1] = t1;
189
 
      cblas_scal(M-2, ATL_rone/t1, A1+2, 1);
190
 
   }
191
 
   else
192
 
   {
193
 
      if (imax != -1) ipiv[1] = imax;
194
 
      else ipiv[1] = 1;
195
 
      if (!iret) iret = 2;
196
 
   }
197
 
   return(iret);
198
 
}
199
 
#endif
200
154
#define MySwap(N_, A_, lda_, ip_) \
201
155
{ \
202
156
   TYPE t0_, t1_, t2_, t3_; \
312
266
 */
313
267
{
314
268
   const int MN = Mmin(M, N);
315
 
   int Nleft, Nright, i, ierr=0;
 
269
   int Nleft, Nright, k, i, ierr=0;
316
270
   #ifdef TCPLX
317
271
      const TYPE one[2] = {ATL_rone, ATL_rzero};
318
272
      const TYPE none[2] = {ATL_rnone, ATL_rzero};
324
278
   #endif
325
279
   TYPE *Ac, *An;
326
280
 
 
281
   if (((size_t)M)*N <= ATL_L1elts)
 
282
      return(Mjoin(PATL,getf2)(M, N, A, lda, ipiv));
 
283
   #if defined(ATL_USEPTHREADS) && defined(ATL_USEPCA)
 
284
      if (N <= (NB<<2) && N >= 16 && M-N >= ATL_PCAMin &&
 
285
          ((size_t)ATL_MulBySize(M)*N) <= CacheEdge*ATL_NTHREADS)
 
286
      {
 
287
         if (N >= 16)
 
288
            ierr = Mjoin(PATL,tgetf2)(M, N, A, lda, ipiv);
 
289
         else
 
290
            ierr = Mjoin(PATL,tgetf2_nocp)(M, N, A, lda, ipiv);
 
291
         return(ierr);
 
292
      }
 
293
   #endif
327
294
   if (MN > ATL_luMmin)
328
295
   {
329
296
      Nleft = MN >> 1;
371
338
         tmp = A[i];
372
339
         if (tmp != ATL_rzero)
373
340
         {
374
 
            cblas_scal(M, ATL_rone/tmp, A, 1);
 
341
            if (Mabs(tmp) > ATL_laSAFMIN)
 
342
               cblas_scal(M, ATL_rone/tmp, A, 1);
 
343
            else
 
344
               for (k=0; k < N; k++)
 
345
                  A[k] /= tmp;
375
346
            A[i] = *A;
376
347
            *A = tmp;
377
348
         }
382
353
         tmp[1] = A[i+1];
383
354
         if (tmp[0] != ATL_rzero || tmp[1] != ATL_rzero)
384
355
         {
385
 
            ATL_CplxInv(tmp, inv);
386
 
            cblas_scal(M, inv, A, 1);
 
356
            if (ATL_lapy2(tmp[0], tmp[1]) >= ATL_laSAFMIN)
 
357
            {
 
358
               ATL_CplxInv(tmp, inv);
 
359
               cblas_scal(M, inv, A, 1);
 
360
            }
 
361
            else
 
362
               Mjoin(PATL,cplxdivide)(M, tmp, A, 1, A, 1);
387
363
            A[i] = *A;
388
364
            A[i+1] = A[1];
389
365
            *A = tmp[0];