← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/lapack/ATL_ormql.c

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot, Sylvestre Ledru, Sébastien Villemot
  • Date: 2013-06-11 15:58:16 UTC
  • mfrom: (1.1.4) (25 sid)
  • mto: This revision was merged to the branch mainline in revision 26.
  • Revision ID: package-import@ubuntu.com-20130611155816-8xeeiziu1iml040c
Tags: 3.10.1-1
http://bugs.debian.org/609287
http://bugs.debian.org/602524
http://bugs.debian.org/701921
http://bugs.debian.org/666203
http://bugs.debian.org/701068
http://bugs.debian.org/697431
[ Sylvestre Ledru ]
* New upstream release (Closes: #609287)

[ Sébastien Villemot ]
* Provide architectural defaults (i.e. precomputed timings) for all
  release archs (except armel and mips for the time being, due to slow
  porterboxes). This will make the package build much faster and should
  eliminate transient build failures due to excessive variance in the
  timings.
* Move symlinks for lib{cblas,f77blas,atlas,lapack_atlas} out of the
  libblas.so.3 alternative and make them always present, so that
  software relying on these libs do not break when another alternative
  is selected for BLAS
* ATLAS now has improved ARM support with native asm constructs. This required
  the following tunes:
  + armel-is-v4t.diff: new patch, prevents FTBFS on armel; otherwise,
    ATLAS uses asm constructs too recent for the platform (armel is only v4t)
  + debian/rules: on armhf, define the ATL_ARM_HARDFP flag; otherwise the asm
    constructs use the soft-float ABI for passing floating points
  + on armhf, ensure that -mfloat-abi=softfp and -mcpu=vfpv3 flags are never
    used; this is implemented via a patch (armhf.diff) and by the use of fixed
    archdefs
* The generic package is now built without multi-threading, because otherwise
  the package fails to build on some single-processor machines (this required
  the introduction of a patch: fix-non-threaded-build.diff). As a side effect,
  the build of the custom package gracefully handles non-threaded
  builds. (Closes: #602524)
* Add libblas.a as slave in the libblas.so alternative (Closes: #701921)
* Add symlinks for lib{f77blas,atlas}.a in /usr/lib (Closes: #666203)
* Modify shlibs file of libatlas3-base, such that packages using
  libblas/liblapack depend on any BLAS/LAPACK alternative, while packages
  depending on ATLAS-specific libraries (e.g. libatlas.so) depend specifically
  on libatlas3-base.
* corei1.diff: remove patch, applied upstream
* Use my @debian.org email address
* Remove obsolete DM-Upload-Allowed flag
* Switch VCS to git
* Remove Conflicts/Replaces against pre-squeeze packages
* libatlas-base-dev now provides libblas.so, as libblas-dev
* No longer use -Wa,--noexecstack in CFLAGS, it makes the package FTBFS
* Do not use POWER3 arch for powerpcspe port (Closes: #701068)
* Bump to debhelper compat level 9
* README.Debian: mention that devscripts is needed to compile the custom
  package (Closes: #697431)
* Bump Standards-Version to 3.9.4. As a consequence, add Built-Using
  fields because the package embeds stuff from liblapack-pic

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lapack/ATL_ormql.c
 
1
 
 
2
#include "stdio.h"
 
3
#include "cblas.h"
 
4
#include "atlas_lapack.h"
 
5
 
 
6
#ifdef  SREAL
 
7
   #define MYOPT LASreal
 
8
#endif
 
9
#ifdef  DREAL
 
10
   #define MYOPT  LADreal
 
11
#endif
 
12
#ifdef  SCPLX
 
13
   #define MYOPT  LAScplx
 
14
#endif
 
15
#ifdef  DCPLX
 
16
   #define MYOPT  LADcplx
 
17
#endif
 
18
 
 
19
int ATL_ormql
 
20
   (const enum CBLAS_SIDE SIDE, const enum CBLAS_TRANSPOSE TRANS,
 
21
    ATL_CINT M, ATL_CINT N, ATL_CINT K, TYPE *A, ATL_CINT lda,
 
22
    const TYPE * TAU, TYPE *C, ATL_CINT ldc, TYPE *WORK, ATL_CINT LWORK)
 
23
/*
 
24
 * This is the C translation of the standard LAPACK Fortran routine:
 
25
 *      SUBROUTINE ATL_ormql( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
 
26
 *                        WORK, LWORK, INFO )
 
27
 *
 
28
 * ATL_ormql.c :
 
29
 * int ATL_ormql(const enum CBLAS_SIDE SIDE SIDE,
 
30
 *        const enum CBLAS_TRANSPOSE TRANS, ATL_CINT M, ATL_CINT N,
 
31
 *        ATL_CINT K, TYPE * A, ATL_CINT lda,TYPE * TAU, TYPE * C, ATL_CINT ldc,
 
32
 *                       TYPE * WORK, ATL_CINT LWORK)
 
33
 *
 
34
 *      NOTE :   ATL_ormql.c will get compiled to four precisions
 
35
 *               single precision real,      double precision real
 
36
 *               single precision complex,   double precision complex
 
37
 *
 
38
 *
 
39
 *  Purpose
 
40
 *  =======
 
41
 *
 
42
 *  ATL_ormql overwrites the general real M-by-N matrix C with
 
43
 *
 
44
 *                  SIDE = 'L'     SIDE = 'R'
 
45
 *  TRANS = 'N':      Q * C          C * Q
 
46
 *  TRANS = 'T':      Q**T * C       C * Q**T
 
47
 *
 
48
 *  where Q is,
 
49
 *        a real orthogonal matrix defined as the product of k
 
50
 *        elementary reflectors
 
51
 *
 
52
 *        Q = H(k) . . . H(2) H(1)
 
53
 *
 
54
 *   OR
 
55
 *        a complex unitary matrix defined as a product of k
 
56
 *        elementary reflectors
 
57
 *
 
58
 *        Q = H(k) . . . H(2) H(1)
 
59
 *
 
60
 *  as returned by ATLL_geqrf.c. Q is of order M if SIDE = 'L' and of order N
 
61
 *  if SIDE = 'R'.
 
62
 *
 
63
 *  Arguments
 
64
 *  =========
 
65
 *
 
66
 *  SIDE    (input) CHARACTER*1
 
67
 *          = 'L': apply Q or Q**T from the Left;
 
68
 *          = 'R': apply Q or Q**T from the Right.
 
69
 *
 
70
 *  TRANS   (input) CHARACTER*1
 
71
 *          = 'N':  No transpose, apply Q;
 
72
 *          = 'T':  Transpose, apply Q**T.
 
73
 *
 
74
 *  M       (input) INTEGER
 
75
 *          The number of rows of the matrix C. M >= 0.
 
76
 *
 
77
 *  N       (input) INTEGER
 
78
 *          The number of columns of the matrix C. N >= 0.
 
79
 *
 
80
 *  K       (input) INTEGER
 
81
 *          The number of elementary reflectors whose product defines
 
82
 *          the matrix Q.
 
83
 *          If SIDE = 'L', M >= K >= 0;
 
84
 *          if SIDE = 'R', N >= K >= 0.
 
85
 *
 
86
 *  A       (input) array, dimension (LDA,K)
 
87
 *          The i-th column must contain the vector which defines the
 
88
 *          elementary reflector H(i), for i = 1,2,...,k, as returned by
 
89
 *          DGEQRF in the first k columns of its array argument A.
 
90
 *          A is modified by the routine but restored on exit.
 
91
 *
 
92
 *  lda     (input) INTEGER
 
93
 *          The leading dimension of the array A.
 
94
 *          If SIDE = 'L', LDA >= max(1,M);
 
95
 *          if SIDE = 'R', LDA >= max(1,N).
 
96
 *
 
97
 *  TAU     (input)  array, dimension (K)
 
98
 *          TAU(i) must contain the scalar factor of the elementary
 
99
 *          reflector H(i), as returned by ATL_geqrf.c.
 
100
 *
 
101
 *  C       (input/output)  array, dimension (LDC,N)
 
102
 *          On entry, the M-by-N matrix C.
 
103
 *          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
 
104
 *
 
105
 *  ldc     (input) INTEGER
 
106
 *          The leading dimension of the array C. LDC >= max(1,M).
 
107
 *
 
108
 *  WORK    (workspace/output) array, dimension (MAX(1,LWORK))
 
109
 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 
110
 *
 
111
 *  LWORK   (input) INTEGER
 
112
 *          The dimension of the array WORK.
 
113
 *          If SIDE = 'L', LWORK >= max(1,N);
 
114
 *          if SIDE = 'R', LWORK >= max(1,M).
 
115
 *          For optimum performance LWORK >= N*NB if SIDE = 'L', and
 
116
 *          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
 
117
 *          blocksize.
 
118
 *
 
119
 *          If LWORK = -1, then a workspace query is assumed; the routine
 
120
 *          only calculates the optimal size of the WORK array, returns
 
121
 *          this value as the first entry of the WORK array, and no error
 
122
 *          message related to LWORK is issued by XERBLA.
 
123
 *
 
124
 *  INFO    (output) INTEGER
 
125
 *          = 0:  successful exit
 
126
 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 
127
 *
 
128
 */
 
129
{
 
130
   ATL_CINT minMN = Mmin(M, N), maxMN = Mmax(M, N);
 
131
   ATL_INT n, nb, j, ib, mi, ni, ic, jc ;
 
132
   TYPE  *ws_QL2,  *ws_T, *ws_larfb;        /* Workspace for QR2,T, larfb     */
 
133
   void *vp=NULL;
 
134
 
 
135
   nb = clapack_ilaenv(LAIS_OPT_NB, LAormqr, MYOPT+LARight+LALower, M, N,K,-1);
 
136
 
 
137
/*
 
138
 * If it is a workspace query, return the size of work required.
 
139
 *    wrksz = wrksz of ATL_larfb + ATL_larft + ATL_geqr2
 
140
 */
 
141
   if (LWORK < 0)
 
142
   {
 
143
      if(SIDE == CblasLeft)
 
144
      {
 
145
         *WORK = ( N*nb + nb*nb + maxMN )  ;
 
146
      }
 
147
      else
 
148
      {
 
149
         *WORK = ( M*nb + nb*nb + maxMN )  ;
 
150
      }
 
151
      return(0);
 
152
   }
 
153
   else if (M < 1 || N < 1)                 /* quick return if no work to do  */
 
154
      return(0);
 
155
/*
 
156
 * If the user gives us too little space, see if we can allocate it ourselves
 
157
 */
 
158
   else
 
159
   {
 
160
      if(SIDE == CblasLeft)
 
161
      {
 
162
         if (LWORK < (N*nb + nb*nb + maxMN))
 
163
         {
 
164
            vp = malloc(ATL_MulBySize(N*nb + nb*nb + maxMN) + ATL_Cachelen);
 
165
            if (!vp)
 
166
               return(-7);
 
167
            WORK = ATL_AlignPtr(vp);
 
168
         }
 
169
      }
 
170
      else
 
171
      {
 
172
         if (LWORK < (M*nb + nb*nb + maxMN))
 
173
         {
 
174
            vp = malloc(ATL_MulBySize(M*nb + nb*nb + maxMN) + ATL_Cachelen);
 
175
            if (!vp)
 
176
               return(-7);
 
177
            WORK = ATL_AlignPtr(vp);
 
178
         }
 
179
      }                                     /* if CblasRight                  */
 
180
   }                                        /* if else                        */
 
181
 
 
182
/*
 
183
 * Assign workspace areas for ATL_larft, ATL_geql2, ATL_larfb
 
184
 */
 
185
 
 
186
   ws_T = WORK;                             /* T at begining of work          */
 
187
   ws_QL2 = WORK +(nb SHIFT)*nb;            /* After T Work space             */
 
188
   ws_larfb = ws_QL2 + (maxMN SHIFT);       /* After workspace for T and QR2  */
 
189
 
 
190
 
 
191
   if (SIDE == CblasLeft)
 
192
   {
 
193
      if (TRANS == CblasNoTrans)
 
194
      {
 
195
         for (j=0 ; j < K ; j = j + nb)
 
196
         {
 
197
            ib = Mmin(K-j, nb);
 
198
/*
 
199
 *          Form the triangular factor of the block reflector
 
200
 *          H = H(i+ib-1) . . . H(i+1) H(i)
 
201
 */
 
202
            ATL_larft(LABackward, LAColumnStore, M-K+j+ib, ib,
 
203
                      A+(j SHIFT)*lda , lda, TAU+(j SHIFT),
 
204
                                   ws_T, ib);
 
205
/*
 
206
 *          H or H' is applied to C(1:m-k+i+ib-1,1:n)
 
207
 */
 
208
            mi = M -K + j +ib ;
 
209
 
 
210
            ATL_larfb(SIDE, TRANS, LABackward, LAColumnStore,
 
211
                      mi, N, ib, A+(j SHIFT)*lda, lda, ws_T, ib,
 
212
                      C, ldc, ws_larfb, N);
 
213
         }
 
214
      }                                     /* CblasNoTrans                   */
 
215
      else
 
216
      {
 
217
         j = (K/nb)*nb;
 
218
         if (j == K)
 
219
         {
 
220
            j=K -nb;
 
221
         }
 
222
 
 
223
         for (j; j >= 0; j = j - nb)
 
224
         {
 
225
            ib = nb;
 
226
            if ((j+nb) > K)
 
227
            {
 
228
               ib = K - j;
 
229
            }
 
230
/*
 
231
 *          Form the triangular factor of the block reflector
 
232
 *          H = H(i+ib-1) . . . H(i+1) H(i)
 
233
 */
 
234
            ATL_larft(LABackward, LAColumnStore, M-K+j+ib, ib,
 
235
                      A+(j SHIFT)*(lda), lda, TAU+(j SHIFT), ws_T, ib);
 
236
/*
 
237
 *          H or H' is applied to C(1:m-k+i+ib-1,1:n)
 
238
 */
 
239
            mi = M -K + j +ib ;
 
240
 
 
241
            ATL_larfb(SIDE, TRANS, LABackward, LAColumnStore,
 
242
                      mi, N, ib, A+(j SHIFT)*(lda), lda, ws_T, ib,
 
243
                      C, ldc, ws_larfb, N);
 
244
         }                                  /* for                            */
 
245
      }
 
246
   }                                        /* cblasLeft main logic           */
 
247
   else
 
248
   {
 
249
      if ( TRANS == CblasNoTrans )
 
250
      {
 
251
         j = (K/nb)*nb;
 
252
         if (j == K)
 
253
         {
 
254
            j=K -nb;
 
255
         }
 
256
         for (j; j >= 0; j = j - nb)
 
257
         {
 
258
 
 
259
            ib = nb;
 
260
            if ((j+nb) > K)
 
261
            {
 
262
               ib = K - j;
 
263
            }
 
264
/*
 
265
 *          Form the triangular factor of the block reflector
 
266
 *          H = H(i+ib-1) . . . H(i+1) H(i)
 
267
 */
 
268
            ATL_larft(LABackward, LAColumnStore, N-K+j+ib, ib,
 
269
                      A+(j SHIFT)*(lda), lda, TAU+(j SHIFT), ws_T, ib);
 
270
/*
 
271
 *          H or H' is applied to C(1:m-k+i+ib-1,1:n)
 
272
 */
 
273
            ni = N - K + j +ib;
 
274
 
 
275
            ATL_larfb(SIDE, TRANS, LABackward, LAColumnStore,
 
276
                      M, ni, ib, A+(j SHIFT)*(lda), lda, ws_T, ib,
 
277
                      C, ldc, ws_larfb, M);
 
278
         }                                  /* for                            */
 
279
 
 
280
      }                                     /* Cblas Tran on Right            */
 
281
      else
 
282
      {
 
283
         for (j = 0 ; j < K; j = j + nb)
 
284
         {
 
285
            ib = Mmin(K-j, nb);
 
286
/*
 
287
 *          Form the triangular factor of the block reflector
 
288
 *          H = H(i+ib-1) . . . H(i+1) H(i)
 
289
 */
 
290
            ATL_larft(LABackward, LAColumnStore, N-K+j+ib, ib,
 
291
                      A+(j SHIFT)*lda , lda, TAU+(j SHIFT),
 
292
                                   ws_T, ib);
 
293
/*
 
294
 *          H or H' is applied to C(1:m-k+i+ib-1,1:n)
 
295
 */
 
296
            ni = N - K + j +ib;
 
297
            ATL_larfb(SIDE, TRANS, LABackward, LAColumnStore,
 
298
                      M, ni, ib, A+(j SHIFT)*lda, lda, ws_T, ib,
 
299
                      C, ldc, ws_larfb, M);
 
300
         }                                  /* for                            */
 
301
      }                                     /* Cblastrans                     */
 
302
   }                                        /* CblasRight                     */
 
303
 
 
304
   if (vp)
 
305
      free(vp);
 
306
   return(0);
 
307
}                                           /* END ATL_ormql                  */
 
308