2
* Automatically Tuned Linear Algebra Software v3.10.1
3
* Copyright (C) 2009 Siju Samuel
5
* Code contributers : Siju Samuel, Anthony M. Castaldo, R. Clint Whaley
7
* Redistribution and use in source and binary forms, with or without
8
* modification, are permitted provided that the following conditions
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.
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.
34
* This is the C translation of the standard LAPACK Fortran routine:
35
* SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO )
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
44
* b) This routine will not validate the input parameters.
48
* ATL_geqr2 computes a QR factorization of a real/complex m by n matrix A:
55
* The number of rows of the matrix A. M >= 0.
58
* The number of columns of the matrix A. N >= 0.
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).
70
* The leading dimension of the array A. LDA >= max(1,M).
72
* TAU (output) array, dimension (min(M,N))
73
* The scalar factors of the elementary reflectors (see Further
76
* WORK (workspace) array, dimension (N)
78
* INFO (output) INTEGER
79
* = 0: successful exit
80
* < 0: if INFO = -i, the i-th argument had an illegal value
85
* The matrix Q is represented as a product of elementary reflectors
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
91
* H(i) = I - tau * v * v' (for Real Precision)
92
* H(i) = I - tau * v * conjugate(v)' (for Complex Precision)
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),
99
#include "atlas_misc.h"
101
#include "atlas_lapack.h"
104
int ATL_geqr2(ATL_CINT M, ATL_CINT N, TYPE *A, ATL_CINT LDA, TYPE *TAU,
107
int lda2 = LDA SHIFT;
110
const TYPE ONE = ATL_rone;
114
const TYPE ONE[2] = {ATL_rone, ATL_rzero};
126
* Generate elementary reflector H(i) to annihilate A(i+1:m,i)
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)) );
132
/* If not last column */
133
if (i < (N-1)) /* If not last column, */
136
* Apply H(i) to A(i:m,i+1:n) from the left
143
AII[0] = A[(i SHIFT)+i*lda2];
144
AII[1] = A[(i SHIFT)+i*lda2 + 1];
146
A[(i SHIFT)+i*lda2] = ONE[0];
147
A[(i SHIFT)+i*lda2 + 1] = ONE[1];
149
TAUVAL[0] = TAU[i SHIFT];
151
TAUVAL[1] = 0.0 -TAU[(i SHIFT) + 1];
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);
157
/* Reassign the values of A[i] */
159
A[(i SHIFT)+i*lda2] = AII;
161
A[(i SHIFT) +i*lda2] = AII[0];
162
A[(i SHIFT) +i*lda2 + 1] = AII[1];
165
} /* END for each column. */
167
} /* END ATL_geqr2 */