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

« back to all changes in this revision

Viewing changes to src/lapack/ATL_geqr2.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:
 
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
/*
 
34
 * This is the C translation of the standard LAPACK Fortran routine:
 
35
 *      SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO )
 
36
 *
 
37
 * ATL_geqr2.c :
 
38
 * int ATL_geqr2( const int M, const int N, TYPE *A, int LDA,
 
39
 *                                                      TYPE  *TAU, TYPE *WORK)
 
40
 *     NOTE :a)   ATL_geqr2.c will get compiled to four precisions
 
41
 *               single precision real,      double precision real
 
42
 *               single precision complex,   double precision complex
 
43
 *
 
44
 *           b) This routine will not validate the input parameters.
 
45
 *  Purpose
 
46
 *  =======
 
47
 *
 
48
 *  ATL_geqr2  computes a QR factorization of a real/complex  m by n matrix A:
 
49
 *  A = Q * R.
 
50
 *
 
51
 *  Arguments
 
52
 *  =========
 
53
 *
 
54
 *  M       (input) INTEGER
 
55
 *          The number of rows of the matrix A.  M >= 0.
 
56
 *
 
57
 *  N       (input) INTEGER
 
58
 *          The number of columns of the matrix A.  N >= 0.
 
59
 *
 
60
 *  A       (input/output) array, dimension (LDA,N)
 
61
 *          On entry, the m by n matrix A.
 
62
 *          On exit, the elements on and above the diagonal of the array
 
63
 *          contain the min(m,n) by n upper trapezoidal matrix R (R is
 
64
 *          upper triangular if m >= n); the elements below the diagonal,
 
65
 *          with the array TAU, represent the orthogonal matrix Q
 
66
 *          (unitary matrix incase of complex precision )  as a
 
67
 *          product of elementary reflectors (see Further Details).
 
68
 *
 
69
 *  LDA     (input) INTEGER
 
70
 *          The leading dimension of the array A.  LDA >= max(1,M).
 
71
 *
 
72
 *  TAU     (output) array, dimension (min(M,N))
 
73
 *          The scalar factors of the elementary reflectors (see Further
 
74
 *          Details).
 
75
 *
 
76
 *  WORK    (workspace)  array, dimension (N)
 
77
 *
 
78
 *  INFO    (output) INTEGER
 
79
 *          = 0: successful exit
 
80
 *          < 0: if INFO = -i, the i-th argument had an illegal value
 
81
 *
 
82
 *  Further Details
 
83
 *  ===============
 
84
 *
 
85
 *  The matrix Q is represented as a product of elementary reflectors
 
86
 *
 
87
 *     Q = H(1) H(2) . . . H(k), where k = min(m,n).
 
88
 *                                             (for Real/Complex Precision)
 
89
 *  Each H(i) has the form
 
90
 *
 
91
 *     H(i) = I - tau * v * v'                 (for Real Precision)
 
92
 *     H(i) = I - tau * v * conjugate(v)'      (for Complex  Precision)
 
93
 *
 
94
 *  where tau is a real/complex  scalar, and v is a real/complex vector with
 
95
 *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
 
96
 *  and tau in TAU(i).
 
97
 */
 
98
 
 
99
#include "atlas_misc.h"
 
100
#include "cblas.h"
 
101
#include "atlas_lapack.h"
 
102
 
 
103
 
 
104
int ATL_geqr2(ATL_CINT M, ATL_CINT N, TYPE *A, ATL_CINT LDA, TYPE  *TAU,
 
105
              TYPE *WORK)
 
106
{
 
107
   int lda2 = LDA  SHIFT;
 
108
 
 
109
   #ifdef TREAL
 
110
      const TYPE ONE = ATL_rone;
 
111
      TYPE AII ;
 
112
      TYPE TAUVAL ;
 
113
   #else
 
114
      const TYPE ONE[2] = {ATL_rone, ATL_rzero};
 
115
      TYPE AII[2];
 
116
      TYPE TAUVAL[2] ;
 
117
   #endif
 
118
 
 
119
   int i, k;
 
120
 
 
121
   k = (M < N)?M:N;
 
122
 
 
123
   for (i=0; i<k; i++)
 
124
   {
 
125
/*
 
126
 *    Generate elementary reflector H(i) to annihilate A(i+1:m,i)
 
127
 */
 
128
      int t=((i+1)<(M-1))?(i+1):(M-1);
 
129
      ATL_larfg((M-i), (A+(i SHIFT)+i*lda2),
 
130
                (A+(t SHIFT)+i*lda2), 1, (TAU+(i SHIFT)) );
 
131
 
 
132
/*    If not last column                   */
 
133
      if (i < (N-1))                        /* If not last column,            */
 
134
      {
 
135
/*
 
136
 *       Apply H(i) to A(i:m,i+1:n) from the left
 
137
 */
 
138
         #ifdef TREAL
 
139
            AII = A[i+i*lda2];
 
140
            A[i+i*lda2] = ONE;
 
141
            TAUVAL = TAU[i];
 
142
         #else
 
143
            AII[0] = A[(i SHIFT)+i*lda2];
 
144
            AII[1] = A[(i SHIFT)+i*lda2 + 1];
 
145
 
 
146
            A[(i SHIFT)+i*lda2] = ONE[0];
 
147
            A[(i SHIFT)+i*lda2 + 1] = ONE[1];
 
148
 
 
149
            TAUVAL[0] = TAU[i SHIFT];
 
150
/*          Conjugate                 */
 
151
            TAUVAL[1] = 0.0 -TAU[(i SHIFT) + 1];
 
152
         #endif
 
153
 
 
154
         ATL_larf(CblasLeft, M-i, N-i-1, (A+(i SHIFT)+i*lda2), 1, TAUVAL ,
 
155
                  (A+(i SHIFT)+(i+1)*lda2) , LDA, WORK);
 
156
 
 
157
         /*  Reassign the values of A[i] */
 
158
         #ifdef TREAL
 
159
            A[(i SHIFT)+i*lda2] = AII;
 
160
         #else
 
161
            A[(i SHIFT) +i*lda2] = AII[0];
 
162
            A[(i SHIFT) +i*lda2 + 1] = AII[1];
 
163
         #endif
 
164
      }
 
165
   } /* END for each column. */
 
166
   return(0);
 
167
} /* END ATL_geqr2 */