← Back to branch summary

~ubuntu-branches/ubuntu/trusty/atlas/trusty

« back to all changes in this revision

Viewing changes to src/lapack/ATL_ormqr.c

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-07-27 14:26:05 UTC
  • mfrom: (18.1.8 sid)
  • Revision ID: package-import@ubuntu.com-20130727142605-5rh3p972h1whdo99
Tags: 3.10.1-2
* Allow the generic package to build on machines with CPU throttling
  enabled. Otherwise the package FTBFS on some buildds (e.g. biber).
  Implementation is done by reactivating the "-Si cputhrchk 0" flag
  (cpu-throtthling-check.diff), and using it in debian/rules.
* Add architectural defaults for armel and mips.
* armhf.diff: do not enforce 32-registers FPU for Fortran

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lapack/ATL_ormqr.c
 
1
/*
 
2
 *             Automatically Tuned Linear Algebra Software v3.10.1
 
3
 * Copyright (C) 2009 Siju Samuel
 
4
 *
 
5
 * Code contributers : Siju Samuel, Anthony M. Castaldo, R. Clint Whaley
 
6
 *
 
7
 * Redistribution and use in source and binary forms, with or without
 
8
 * modification, are permitted provided that the following conditions
 
9
 * are met:
 
10
 *   1. Redistributions of source code must retain the above copyright
 
11
 *      notice, this list of conditions and the following disclaimer.
 
12
 *   2. Redistributions in binary form must reproduce the above copyright
 
13
 *      notice, this list of conditions, and the following disclaimer in the
 
14
 *      documentation and/or other materials provided with the distribution.
 
15
 *   3. The name of the ATLAS group or the names of its contributers may
 
16
 *      not be used to endorse or promote products derived from this
 
17
 *      software without specific written permission.
 
18
 *
 
19
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 
20
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 
21
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 
22
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
 
23
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 
24
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 
25
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 
26
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 
27
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 
28
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 
29
 * POSSIBILITY OF SUCH DAMAGE.
 
30
 *
 
31
 */
 
32
 
 
33
#include "stdio.h"
 
34
#include "cblas.h"
 
35
#include "atlas_lapack.h"
 
36
 
 
37
#ifdef  SREAL
 
38
     #define MYOPT LASreal
 
39
#endif
 
40
#ifdef  DREAL
 
41
    #define MYOPT  LADreal
 
42
#endif
 
43
#ifdef  SCPLX
 
44
    #define MYOPT  LAScplx
 
45
#endif
 
46
#ifdef  DCPLX
 
47
    #define MYOPT  LADcplx
 
48
#endif
 
49
 
 
50
int ATL_ormqr
 
51
   (const enum CBLAS_SIDE SIDE, const enum CBLAS_TRANSPOSE TRANS,
 
52
    ATL_CINT M, ATL_CINT N, ATL_CINT K, TYPE *A, ATL_CINT lda,
 
53
    const TYPE *TAU, TYPE *C, ATL_CINT ldc, TYPE *WORK, ATL_CINT LWORK)
 
54
/*
 
55
 * This is the C translation of the standard LAPACK Fortran routine:
 
56
 *      SUBROUTINE ATL_ormqr( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
 
57
 *                        WORK, LWORK, INFO )
 
58
 *
 
59
 * ATL_ormqr.c :
 
60
 * int ATL_ormqr(const enum CBLAS_SIDE SIDE SIDE,
 
61
 *        const enum CBLAS_TRANSPOSE TRANS, ATL_CINT M, ATL_CINT N,
 
62
 *        ATL_CINT K, TYPE * A, ATL_CINT lda,TYPE * TAU, TYPE * C, ATL_CINT ldc,
 
63
 *                       TYPE * WORK, ATL_CINT LWORK)
 
64
 *
 
65
 *      NOTE :   ATL_ormqr.c will get compiled to four precisions
 
66
 *               single precision real,      double precision real
 
67
 *               single precision complex,   double precision complex
 
68
 *
 
69
 *
 
70
 *
 
71
 *  Purpose
 
72
 *  =======
 
73
 *
 
74
 *  ATL_ormqr overwrites the general real M-by-N matrix C with
 
75
 *
 
76
 *                  SIDE = 'L'     SIDE = 'R'
 
77
 *  TRANS = 'N':      Q * C          C * Q
 
78
 *  TRANS = 'T':      Q**T * C       C * Q**T
 
79
 *
 
80
 *  where Q is,
 
81
 *        a real orthogonal matrix defined as the product of k
 
82
 *        elementary reflectors
 
83
 *
 
84
 *        Q = H(1) H(2) . . . H(k)
 
85
 *
 
86
 *   OR
 
87
 *        a complex unitary matrix defined as a product of k
 
88
 *        elementary reflectors
 
89
 *
 
90
 *        Q = H(1) H(2) . . . H(k)
 
91
 *
 
92
 *  as returned by ATLL_geqrf.c. Q is of order M if SIDE = 'L' and of order N
 
93
 *  if SIDE = 'R'.
 
94
 *
 
95
 *  Arguments
 
96
 *  =========
 
97
 *
 
98
 *  SIDE    (input) CHARACTER*1
 
99
 *          = 'L': apply Q or Q**T from the Left;
 
100
 *          = 'R': apply Q or Q**T from the Right.
 
101
 *
 
102
 *  TRANS   (input) CHARACTER*1
 
103
 *          = 'N':  No transpose, apply Q;
 
104
 *          = 'T':  Transpose, apply Q**T.
 
105
 *
 
106
 *  M       (input) INTEGER
 
107
 *          The number of rows of the matrix C. M >= 0.
 
108
 *
 
109
 *  N       (input) INTEGER
 
110
 *          The number of columns of the matrix C. N >= 0.
 
111
 *
 
112
 *  K       (input) INTEGER
 
113
 *          The number of elementary reflectors whose product defines
 
114
 *          the matrix Q.
 
115
 *          If SIDE = 'L', M >= K >= 0;
 
116
 *          if SIDE = 'R', N >= K >= 0.
 
117
 *
 
118
 *  A       (input) array, dimension (LDA,K)
 
119
 *          The i-th column must contain the vector which defines the
 
120
 *          elementary reflector H(i), for i = 1,2,...,k, as returned by
 
121
 *          DGEQRF in the first k columns of its array argument A.
 
122
 *          A is modified by the routine but restored on exit.
 
123
 *
 
124
 *  lda     (input) INTEGER
 
125
 *          The leading dimension of the array A.
 
126
 *          If SIDE = 'L', LDA >= max(1,M);
 
127
 *          if SIDE = 'R', LDA >= max(1,N).
 
128
 *
 
129
 *  TAU     (input)  array, dimension (K)
 
130
 *          TAU(i) must contain the scalar factor of the elementary
 
131
 *          reflector H(i), as returned by ATL_geqrf.c.
 
132
 *
 
133
 *  C       (input/output)  array, dimension (LDC,N)
 
134
 *          On entry, the M-by-N matrix C.
 
135
 *          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
 
136
 *
 
137
 *  ldc     (input) INTEGER
 
138
 *          The leading dimension of the array C. LDC >= max(1,M).
 
139
 *
 
140
 *  WORK    (workspace/output) array, dimension (MAX(1,LWORK))
 
141
 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 
142
 *
 
143
 *  LWORK   (input) INTEGER
 
144
 *          The dimension of the array WORK.
 
145
 *          If SIDE = 'L', LWORK >= max(1,N);
 
146
 *          if SIDE = 'R', LWORK >= max(1,M).
 
147
 *          For optimum performance LWORK >= N*NB if SIDE = 'L', and
 
148
 *          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
 
149
 *          blocksize.
 
150
 *
 
151
 *          If LWORK = -1, then a workspace query is assumed; the routine
 
152
 *          only calculates the optimal size of the WORK array, returns
 
153
 *          this value as the first entry of the WORK array, and no error
 
154
 *          message related to LWORK is issued by XERBLA.
 
155
 *
 
156
 *  INFO    (output) INTEGER
 
157
 *          = 0:  successful exit
 
158
 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 
159
 */
 
160
{
 
161
   ATL_CINT minMN = Mmin(M, N), maxMN = Mmax(M, N);
 
162
   ATL_INT n, nb, j, ib, mi, ni, ic, jc ;
 
163
   TYPE  *ws_QR2,  *ws_T, *ws_larfb;        /* Workspace for QR2,T, larfb     */
 
164
   void *vp=NULL;
 
165
 
 
166
   nb = clapack_ilaenv(LAIS_OPT_NB, LAormqr, MYOPT+LARight+LAUpper, M, N, K,-1);
 
167
 
 
168
/*
 
169
 * If it is a workspace query, return the size of work required.
 
170
 *    wrksz = wrksz of ATL_larfb + ATL_larft + ATL_geqr2
 
171
 */
 
172
   if (LWORK < 0)
 
173
   {
 
174
      if(SIDE == CblasLeft)
 
175
      {
 
176
         *WORK = ( N*nb + nb*nb + maxMN )  ;
 
177
      }
 
178
      else
 
179
      {
 
180
         *WORK = ( M*nb + nb*nb + maxMN )  ;
 
181
      }
 
182
      return(0);
 
183
   }
 
184
   else if (M < 1 || N < 1)                 /* quick return if no work to do  */
 
185
      return(0);
 
186
/*
 
187
 * If the user gives us too little space, see if we can allocate it ourselves
 
188
 */
 
189
   else
 
190
   {
 
191
      if(SIDE == CblasLeft)
 
192
      {
 
193
         if (LWORK < (N*nb + nb*nb + maxMN))
 
194
         {
 
195
            vp = malloc(ATL_MulBySize(N*nb + nb*nb + maxMN) + ATL_Cachelen);
 
196
            if (!vp)
 
197
               return(-7);
 
198
            WORK = ATL_AlignPtr(vp);
 
199
         }
 
200
      }
 
201
      else
 
202
      {
 
203
         if (LWORK < (M*nb + nb*nb + maxMN))
 
204
         {
 
205
            vp = malloc(ATL_MulBySize(M*nb + nb*nb + maxMN) + ATL_Cachelen);
 
206
            if (!vp)
 
207
               return(-7);
 
208
            WORK = ATL_AlignPtr(vp);
 
209
         }
 
210
      } /* if CblasRight */
 
211
   }
 
212
 
 
213
/*
 
214
 * Assign workspace areas for ATL_larft, ATL_geqr2, ATL_larfb
 
215
 */
 
216
 
 
217
   ws_T = WORK;                             /* T at begining of work          */
 
218
   ws_QR2 = WORK +(nb SHIFT)*nb;            /* After T Work space             */
 
219
   ws_larfb = ws_QR2 + (maxMN SHIFT);       /* After workspace for T and QR2  */
 
220
 
 
221
 
 
222
   if (SIDE == CblasLeft)
 
223
   {
 
224
      if ( TRANS == CblasNoTrans )
 
225
      {
 
226
         j = (K/nb)*nb;
 
227
         if (j == K)
 
228
         {
 
229
            j=K -nb;
 
230
         }
 
231
         for (j; j >= 0; j = j - nb)
 
232
         {
 
233
            ib = nb;
 
234
            if ((j+nb) > K)
 
235
            {
 
236
               ib = K - j;
 
237
            }
 
238
/*
 
239
 *          Form the triangular factor of the block reflector
 
240
 *          H = H(i) H(i+1) . . . H(i+ib-1)
 
241
 */
 
242
            ATL_larft(LAForward, LAColumnStore, M-j, ib, A+(j SHIFT)*(lda+1),
 
243
                      lda, TAU+(j SHIFT), ws_T, ib);
 
244
/*
 
245
 *          H or H' is applied to C(i:m,1:n)
 
246
 */
 
247
            ATL_larfb(SIDE, TRANS, LAForward, LAColumnStore,
 
248
                      (M-j), N, ib, A+(j SHIFT)*(lda+1), lda, ws_T, ib,
 
249
                      C+(j SHIFT), ldc, ws_larfb, N);
 
250
          }                                 /* for                            */
 
251
      }                                     /* CblasNoTrans                   */
 
252
      else
 
253
      {
 
254
         for (j = 0 ; j < K; j = j + nb)
 
255
         {
 
256
            ib = Mmin(K-j, nb);
 
257
/*
 
258
 *          Form the triangular factor of the block reflector
 
259
 *          H = H(i) H(i+1) . . . H(i+ib-1)
 
260
 */
 
261
            ATL_larft(LAForward, LAColumnStore, M-j, ib, A+(j SHIFT)*(lda+1),
 
262
                      lda, TAU+(j SHIFT), ws_T, ib);
 
263
/*
 
264
 *          H or H' is applied to C(i:m,1:n)
 
265
 */
 
266
            ATL_larfb(SIDE, TRANS, LAForward, LAColumnStore,
 
267
                      (M-j), N, ib, A+(j SHIFT)*(lda+1), lda, ws_T, ib,
 
268
                      C+(j SHIFT), ldc, ws_larfb, N);
 
269
         }                                  /* for                            */
 
270
      }                                     /* CblasNoTran                    */
 
271
   }                                        /* cblasLeft                      */
 
272
   else
 
273
   {
 
274
      if ( TRANS == CblasNoTrans )
 
275
      {
 
276
         for (j = 0 ; j < K; j = j + nb)
 
277
         {
 
278
            ib = Mmin(K-j, nb);
 
279
/*
 
280
 *          Form the triangular factor of the block reflector
 
281
 *          H = H(i) H(i+1) . . . H(i+ib-1)
 
282
 */
 
283
            ATL_larft(LAForward, LAColumnStore, N-j, ib, A+(j SHIFT)*(lda+1),
 
284
                      lda, TAU+(j SHIFT), ws_T, ib);
 
285
/*
 
286
 *              H or H' is applied to C(1:m,i:n)
 
287
 */
 
288
            ATL_larfb(SIDE, TRANS, LAForward, LAColumnStore,
 
289
                      M, N-j, ib, A+(j SHIFT)*(lda+1), lda, ws_T, ib,
 
290
                      C+((j SHIFT)*ldc), ldc, ws_larfb, M);
 
291
          }                                 /* for                            */
 
292
      }
 
293
      else
 
294
      {
 
295
         j = (K/nb)*nb;
 
296
         if (j == K)
 
297
         {
 
298
            j=K -nb;
 
299
         }
 
300
         for (j; j >= 0; j = j - nb)
 
301
         {
 
302
 
 
303
            ib = nb;
 
304
            if ((j+nb) > K)
 
305
            {
 
306
               ib = K - j;
 
307
            }
 
308
/*
 
309
 *          Form the triangular factor of the block reflector
 
310
 *          H = H(i) H(i+1) . . . H(i+ib-1)
 
311
 */
 
312
            ATL_larft(LAForward, LAColumnStore, N-j, ib, A+(j SHIFT)*(lda+1),
 
313
                      lda, TAU+(j SHIFT), ws_T, ib);
 
314
/*
 
315
 *              H or H' is applied to C(1:m,i:n)
 
316
 */
 
317
 
 
318
            ATL_larfb(SIDE, TRANS, LAForward, LAColumnStore,
 
319
                      M, N-j , ib, A+(j SHIFT)*(lda+1), lda, ws_T, ib,
 
320
                      C+((j SHIFT)*ldc) , ldc, ws_larfb, M);
 
321
         }                                  /* for                            */
 
322
 
 
323
      }                                     /* Cblas Tran on Right            */
 
324
   }
 
325
 
 
326
   if (vp)
 
327
      free(vp);
 
328
   return(0);
 
329
}                                           /* END ATL_ormqr                  */
 
330
 
 
331