1
/* ---------------------------------------------------------------------
3
* -- Automatically Tuned Linear Algebra Software (ATLAS)
4
* (C) Copyright 2000 All Rights Reserved
6
* -- ATLAS routine -- Version 3.2 -- December 15, 2000
8
* -- Suggestions, comments, bugs reports should be sent to the follo-
9
* wing e-mail address: atlas@cs.utk.edu
11
* Author : Antoine P. Petitet
12
* University of Tennessee - Innovative Computing Laboratory
13
* Knoxville TN, 37996-1301, USA.
15
* ---------------------------------------------------------------------
17
* -- Copyright notice and Licensing terms:
19
* Redistribution and use in source and binary forms, with or without
20
* modification, are permitted provided that the following conditions
23
* 1. Redistributions of source code must retain the above copyright
24
* notice, this list of conditions and the following disclaimer.
25
* 2. Redistributions in binary form must reproduce the above copyright
26
* notice, this list of conditions, and the following disclaimer in
27
* the documentation and/or other materials provided with the distri-
29
* 3. The name of the University, the ATLAS group, or the names of its
30
* contributors may not be used to endorse or promote products deri-
31
* ved from this software without specific written permission.
35
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
36
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
37
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
38
* A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
39
* OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
40
* CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
41
* TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
42
* OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
43
* RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
44
* CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
45
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
47
* ---------------------------------------------------------------------
52
#include "atlas_refmisc.h"
53
#include "atlas_reflvl2.h"
54
#include "atlas_reflevel2.h"
58
const enum ATLAS_UPLO UPLO,
59
const enum ATLAS_TRANS TRANS,
60
const enum ATLAS_DIAG DIAG,
73
* ATL_creftbsv solves one of the systems of equations
75
* A * x = b, or conjg( A ) * x = b, or
77
* A'* x = b, or conjg( A' ) * x = b,
79
* where b and x are n-element vectors and A is an n by n unit, or non-
80
* unit, upper or lower triangular band matrix, with (k+1) diagonals.
82
* No test for singularity or near-singularity is included in this
83
* routine. Such tests must be performed before calling this routine.
88
* UPLO (input) const enum ATLAS_UPLO
89
* On entry, UPLO specifies whether the matrix is an upper or
90
* lower triangular matrix as follows:
92
* UPLO = AtlasUpper A is an upper triangular matrix.
94
* UPLO = AtlasLower A is a lower triangular matrix.
98
* TRANS (input) const enum ATLAS_TRANS
99
* On entry, TRANS specifies the equations to be solved as fol-
102
* TRANS = AtlasNoTrans A * x = b,
104
* TRANS = AtlasConj conjg( A ) * x = b,
106
* TRANS = AtlasTrans A' * x = b,
108
* TRANS = AtlasConjTrans conjg( A' ) * x = b.
112
* DIAG (input) const enum ATLAS_DIAG
113
* On entry, DIAG specifies whether or not A is unit triangu-
116
* DIAG = AtlasUnit A is assumed to be unit triangular,
118
* DIAG = AtlasNonUnit A is not assumed to be unit trian-
123
* N (input) const int
124
* On entry, N specifies the order of the matrix A. N must be at
125
* least zero. Unchanged on exit.
127
* K (input) const int
128
* On entry with UPLO = AtlasUpper, K specifies the number of
129
* super-diagonals of the matrix A. With UPLO = AtlasLower, K
130
* specifies the number of sub-diagonals of the matrix A. K must
131
* satisfy 0 <= K. Unchanged on exit.
133
* A (input) const float *
134
* On entry, A points to an array of size equal to or greater
135
* than LDA * n * sizeof( float [2] ). Before entry with
136
* UPLO = AtlasUpper, the leading (k + 1) by n part of the array
137
* A must contain the upper triangular band part of the matrix
138
* of coefficients, supplied column by column, with the leading
139
* diagonal of the matrix in row k of the array, the first su-
140
* per-diagonal starting at position 1 in row k-1, and so on.
141
* The top left k by k triangle of the array A is not referen-
142
* ced. The following program segment will transfer an upper
143
* triangular band matrix from conventional full matrix storage
146
* for( j = 0; j < n; j++ )
149
* for( i = ( m < 0 ? -m : 0 ); i < j; i++ )
151
* a[((m+i+j*LDA)<<1)+0] = real( matrix( i, j ) );
152
* a[((m+i+j*LDA)<<1)+1] = imag( matrix( i, j ) );
156
* Before entry with UPLO = AtlasLower, the leading (k + 1) by n
157
* part of the array A must contain the lower triangular band
158
* part of the matrix of coefficients, supplied column by co-
159
* lumn, with the leading diagonal of the matrix in row 0 of the
160
* array, the first sub-diagonal starting at position 0 in row
161
* 1, and so on. The bottom right k by k triangle of the array A
162
* is not referenced. The following program segment will trans-
163
* fer a lower real triangular band matrix from conventional
164
* full matrix storage to band storage:
166
* for( j = 0; j < n; j++ )
168
* i1 = ( n > j + k + 1 ? j + k + 1 : n );
169
* for( i = j; i < i1; i++ )
171
* a[((i-j+j*LDA)<<1)+0] = real( matrix( i, j ) );
172
* a[((i-j+j*LDA)<<1)+1] = imag( matrix( i, j ) );
176
* Note that when DIAG = AtlasUnit the elements of the array A
177
* corresponding to the diagonal elements of the matrix are not
178
* referenced, but are assumed to be unity. Unchanged on exit.
180
* LDA (input) const int
181
* On entry, LDA specifies the leading dimension of A as decla-
182
* red in the calling (sub) program. LDA must be at least
183
* k + 1. Unchanged on exit.
185
* X (input/output) float *
186
* On entry, X points to the first entry to be accessed of an
187
* incremented array of size equal to or greater than
188
* ( 1 + ( n - 1 ) * abs( INCX ) ) * sizeof( float [2] ),
189
* that contains the vector x. Before entry, the incremented ar-
190
* ray X must contain the n element right-hand side vector b. On
191
* exit, X is overwritten with the solution vector x.
193
* INCX (input) const int
194
* On entry, INCX specifies the increment for the elements of X.
195
* INCX must not be zero. Unchanged on exit.
197
* ---------------------------------------------------------------------
200
* .. Executable Statements ..
205
if( UPLO == AtlasUpper )
207
if( TRANS == AtlasNoTrans )
209
if( DIAG == AtlasNonUnit )
211
ATL_creftbsvUNN( N, K, A, LDA, X, INCX );
215
ATL_creftbsvUNU( N, K, A, LDA, X, INCX );
218
else if( TRANS == AtlasConj )
220
if( DIAG == AtlasNonUnit )
222
ATL_creftbsvUCN( N, K, A, LDA, X, INCX );
226
ATL_creftbsvUCU( N, K, A, LDA, X, INCX );
229
else if( TRANS == AtlasTrans )
231
if( DIAG == AtlasNonUnit )
233
ATL_creftbsvUTN( N, K, A, LDA, X, INCX );
237
ATL_creftbsvUTU( N, K, A, LDA, X, INCX );
242
if( DIAG == AtlasNonUnit )
244
ATL_creftbsvUHN( N, K, A, LDA, X, INCX );
248
ATL_creftbsvUHU( N, K, A, LDA, X, INCX );
254
if( TRANS == AtlasNoTrans )
256
if( DIAG == AtlasNonUnit )
258
ATL_creftbsvLNN( N, K, A, LDA, X, INCX );
262
ATL_creftbsvLNU( N, K, A, LDA, X, INCX );
265
else if( TRANS == AtlasConj )
267
if( DIAG == AtlasNonUnit )
269
ATL_creftbsvLCN( N, K, A, LDA, X, INCX );
273
ATL_creftbsvLCU( N, K, A, LDA, X, INCX );
276
else if( TRANS == AtlasTrans )
278
if( DIAG == AtlasNonUnit )
280
ATL_creftbsvLTN( N, K, A, LDA, X, INCX );
284
ATL_creftbsvLTU( N, K, A, LDA, X, INCX );
289
if( DIAG == AtlasNonUnit )
291
ATL_creftbsvLHN( N, K, A, LDA, X, INCX );
295
ATL_creftbsvLHU( N, K, A, LDA, X, INCX );
300
* End of ATL_creftbsv