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_reflvl3.h"
54
#include "atlas_reflevel3.h"
58
const enum ATLAS_UPLO UPLO,
59
const enum ATLAS_TRANS TRANS,
74
* ATL_crefherk performs one of the Hermitian rank k operations
76
* C := alpha * A * conjg( A' ) + beta * C,
80
* C := alpha * conjg( A' ) * A + beta * C,
82
* where alpha and beta are real scalars, C is an n by n Hermitian ma-
83
* trix and A is an n by k matrix in the first case and a k by n matrix
89
* UPLO (input) const enum ATLAS_UPLO
90
* On entry, UPLO specifies whether the upper or lower triangu-
91
* lar part of the array C is to be referenced as follows:
93
* UPLO = AtlasUpper Only the upper triangular part of C
94
* is to be referenced.
96
* UPLO = AtlasLower Only the lower triangular part of C
97
* is to be referenced.
101
* TRANS (input) const enum ATLAS_TRANS
102
* On entry, TRANS specifies the operation to be performed as
105
* TRANS = AtlasNoTrans C := alpha * A*conjg( A )' +
108
* TRANS = AtlasConjTrans C := alpha * conjg( A )'*A +
113
* N (input) const int
114
* On entry, N specifies the order of the matrix C. N must be at
115
* least zero. Unchanged on exit.
117
* K (input) const int
118
* On entry, with TRANS = AtlasNoTrans, K specifies the number
119
* of columns of the matrix A, and otherwise, K specifies the
120
* number of rows of the matrix A. K must be at least zero. Un-
123
* ALPHA (input) const float
124
* On entry, ALPHA specifies the real scalar alpha. When ALPHA
125
* is supplied as zero then the entries of the matrix A need
126
* not be set on input. Unchanged on exit.
128
* A (input) const float *
129
* On entry, A points to an array of size equal to or greater
130
* than LDA * ka * sizeof( float [2] ), where ka is k when
131
* TRANS = AtlasNoTrans, and is n otherwise. Before entry with
132
* TRANS = AtlasNoTrans, the leading n by k part of the array A
133
* must contain the matrix A, otherwise the leading k by n part
134
* of the array A must contain the matrix A. Unchanged on exit.
136
* LDA (input) const int
137
* On entry, LDA specifies the leading dimension of A as decla-
138
* red in the calling (sub) program. LDA must be at least
139
* MAX( 1, n ) when TRANS = AtlasNoTrans, and MAX( 1, k ) other-
140
* wise. Unchanged on exit.
142
* BETA (input) const float
143
* On entry, BETA specifies the real scalar beta. When BETA is
144
* supplied as zero then the entries of the matrix C need not
145
* be set on input. Unchanged on exit.
147
* C (input/output) float *
148
* On entry, C points to an array of size equal to or greater
149
* than LDC * n * sizeof( float [2] ), Before entry with
150
* UPLO = AtlasUpper, the leading n by n upper triangular part
151
* of the array C must contain the upper triangular part of the
152
* Hermitian matrix and the strictly lower triangular part of C
153
* is not referenced. On exit, the upper triangular part of the
154
* array C is overwritten by the upper triangular part of the
155
* updated matrix. Before entry with UPLO = AtlasLower, the
156
* leading n by n lower triangular part of the array C must con-
157
* tain the lower triangular part of the Hermitian matrix and
158
* the strictly upper triangular part of C is not referenced. On
159
* exit, the lower triangular part of the array C is overwritten
160
* by the lower triangular part of the updated matrix.
161
* Note that the imaginary parts of the diagonal elements of C
162
* need not be set, they are assumed to be zero, and on exit
163
* they are set to zero.
165
* LDC (input) const int
166
* On entry, LDC specifies the leading dimension of A as decla-
167
* red in the calling (sub) program. LDC must be at least
168
* MAX( 1, n ). Unchanged on exit.
170
* ---------------------------------------------------------------------
173
* .. Local Variables ..
175
int i, icij, j, jcj, ldc2 = ( LDC << 1 ),
176
ldcp12 = ( ( LDC + 1 ) << 1 );
178
* .. Executable Statements ..
181
if( ( N == 0 ) || ( ( ( ALPHA == ATL_sZERO ) || ( K == 0 ) ) &&
182
( BETA == ATL_sONE ) ) )
185
if( ALPHA == ATL_sZERO )
187
if( UPLO == AtlasUpper )
189
if( BETA == ATL_sZERO )
191
for( j = 0, jcj = 0; j < N; j++, jcj += ldc2 )
193
for( i = 0, icij = jcj; i <= j; i++, icij += 2 )
194
{ Mset( ATL_sZERO, ATL_sZERO, C[icij], C[icij+1] ); }
197
else if( BETA != ATL_sONE )
199
for( j = 0, jcj = 0; j < N; j++, jcj += ldc2 )
201
for( i = 0, icij = jcj; i < j; i++, icij += 2 )
202
{ Mset( BETA * C[icij], BETA * C[icij+1], C[icij], C[icij+1] ); }
203
Mset( BETA * C[icij], ATL_sZERO, C[icij], C[icij+1] );
209
if( BETA == ATL_sZERO )
211
for( j = 0, jcj = 0; j < N; j++, jcj += ldcp12 )
213
for( i = j, icij = jcj; i < N; i++, icij += 2 )
214
{ Mset( ATL_sZERO, ATL_sZERO, C[icij], C[icij+1] ); }
217
else if( BETA != ATL_sONE )
219
for( j = 0, jcj = 0; j < N; j++, jcj += ldcp12 )
221
Mset( BETA * C[jcj], ATL_sZERO, C[jcj], C[jcj+1] );
222
for( i = j+1, icij = jcj+2; i < N; i++, icij += 2 )
223
{ Mset( BETA * C[icij], BETA * C[icij+1], C[icij], C[icij+1] ); }
230
if( UPLO == AtlasUpper )
232
if( TRANS == AtlasNoTrans )
233
{ ATL_crefherkUN( N, K, ALPHA, A, LDA, BETA, C, LDC ); }
235
{ ATL_crefherkUC( N, K, ALPHA, A, LDA, BETA, C, LDC ); }
239
if( TRANS == AtlasNoTrans )
240
{ ATL_crefherkLN( N, K, ALPHA, A, LDA, BETA, C, LDC ); }
242
{ ATL_crefherkLC( N, K, ALPHA, A, LDA, BETA, C, LDC ); }
245
* End of ATL_crefherk