![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLATSP 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE ZLATSP( UPLO, N, X, ISEED ) 00012 * 00013 * .. Scalar Arguments .. 00014 * CHARACTER UPLO 00015 * INTEGER N 00016 * .. 00017 * .. Array Arguments .. 00018 * INTEGER ISEED( * ) 00019 * COMPLEX*16 X( * ) 00020 * .. 00021 * 00022 * 00023 *> \par Purpose: 00024 * ============= 00025 *> 00026 *> \verbatim 00027 *> 00028 *> ZLATSP generates a special test matrix for the complex symmetric 00029 *> (indefinite) factorization for packed matrices. The pivot blocks of 00030 *> the generated matrix will be in the following order: 00031 *> 2x2 pivot block, non diagonalizable 00032 *> 1x1 pivot block 00033 *> 2x2 pivot block, diagonalizable 00034 *> (cycle repeats) 00035 *> A row interchange is required for each non-diagonalizable 2x2 block. 00036 *> \endverbatim 00037 * 00038 * Arguments: 00039 * ========== 00040 * 00041 *> \param[in] UPLO 00042 *> \verbatim 00043 *> UPLO is CHARACTER 00044 *> Specifies whether the generated matrix is to be upper or 00045 *> lower triangular. 00046 *> = 'U': Upper triangular 00047 *> = 'L': Lower triangular 00048 *> \endverbatim 00049 *> 00050 *> \param[in] N 00051 *> \verbatim 00052 *> N is INTEGER 00053 *> The dimension of the matrix to be generated. 00054 *> \endverbatim 00055 *> 00056 *> \param[out] X 00057 *> \verbatim 00058 *> X is COMPLEX*16 array, dimension (N*(N+1)/2) 00059 *> The generated matrix in packed storage format. The matrix 00060 *> consists of 3x3 and 2x2 diagonal blocks which result in the 00061 *> pivot sequence given above. The matrix outside these 00062 *> diagonal blocks is zero. 00063 *> \endverbatim 00064 *> 00065 *> \param[in,out] ISEED 00066 *> \verbatim 00067 *> ISEED is INTEGER array, dimension (4) 00068 *> On entry, the seed for the random number generator. The last 00069 *> of the four integers must be odd. (modified on exit) 00070 *> \endverbatim 00071 * 00072 * Authors: 00073 * ======== 00074 * 00075 *> \author Univ. of Tennessee 00076 *> \author Univ. of California Berkeley 00077 *> \author Univ. of Colorado Denver 00078 *> \author NAG Ltd. 00079 * 00080 *> \date November 2011 00081 * 00082 *> \ingroup complex16_lin 00083 * 00084 * ===================================================================== 00085 SUBROUTINE ZLATSP( UPLO, N, X, ISEED ) 00086 * 00087 * -- LAPACK test routine (version 3.4.0) -- 00088 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00089 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00090 * November 2011 00091 * 00092 * .. Scalar Arguments .. 00093 CHARACTER UPLO 00094 INTEGER N 00095 * .. 00096 * .. Array Arguments .. 00097 INTEGER ISEED( * ) 00098 COMPLEX*16 X( * ) 00099 * .. 00100 * 00101 * ===================================================================== 00102 * 00103 * .. Parameters .. 00104 COMPLEX*16 EYE 00105 PARAMETER ( EYE = ( 0.0D0, 1.0D0 ) ) 00106 * .. 00107 * .. Local Scalars .. 00108 INTEGER J, JJ, N5 00109 DOUBLE PRECISION ALPHA, ALPHA3, BETA 00110 COMPLEX*16 A, B, C, R 00111 * .. 00112 * .. External Functions .. 00113 COMPLEX*16 ZLARND 00114 EXTERNAL ZLARND 00115 * .. 00116 * .. Intrinsic Functions .. 00117 INTRINSIC ABS, SQRT 00118 * .. 00119 * .. Executable Statements .. 00120 * 00121 * Initialize constants 00122 * 00123 ALPHA = ( 1.D0+SQRT( 17.D0 ) ) / 8.D0 00124 BETA = ALPHA - 1.D0 / 1000.D0 00125 ALPHA3 = ALPHA*ALPHA*ALPHA 00126 * 00127 * Fill the matrix with zeros. 00128 * 00129 DO 10 J = 1, N*( N+1 ) / 2 00130 X( J ) = 0.0D0 00131 10 CONTINUE 00132 * 00133 * UPLO = 'U': Upper triangular storage 00134 * 00135 IF( UPLO.EQ.'U' ) THEN 00136 N5 = N / 5 00137 N5 = N - 5*N5 + 1 00138 * 00139 JJ = N*( N+1 ) / 2 00140 DO 20 J = N, N5, -5 00141 A = ALPHA3*ZLARND( 5, ISEED ) 00142 B = ZLARND( 5, ISEED ) / ALPHA 00143 C = A - 2.D0*B*EYE 00144 R = C / BETA 00145 X( JJ ) = A 00146 X( JJ-2 ) = B 00147 JJ = JJ - J 00148 X( JJ ) = ZLARND( 2, ISEED ) 00149 X( JJ-1 ) = R 00150 JJ = JJ - ( J-1 ) 00151 X( JJ ) = C 00152 JJ = JJ - ( J-2 ) 00153 X( JJ ) = ZLARND( 2, ISEED ) 00154 JJ = JJ - ( J-3 ) 00155 X( JJ ) = ZLARND( 2, ISEED ) 00156 IF( ABS( X( JJ+( J-3 ) ) ).GT.ABS( X( JJ ) ) ) THEN 00157 X( JJ+( J-4 ) ) = 2.0D0*X( JJ+( J-3 ) ) 00158 ELSE 00159 X( JJ+( J-4 ) ) = 2.0D0*X( JJ ) 00160 END IF 00161 JJ = JJ - ( J-4 ) 00162 20 CONTINUE 00163 * 00164 * Clean-up for N not a multiple of 5. 00165 * 00166 J = N5 - 1 00167 IF( J.GT.2 ) THEN 00168 A = ALPHA3*ZLARND( 5, ISEED ) 00169 B = ZLARND( 5, ISEED ) / ALPHA 00170 C = A - 2.D0*B*EYE 00171 R = C / BETA 00172 X( JJ ) = A 00173 X( JJ-2 ) = B 00174 JJ = JJ - J 00175 X( JJ ) = ZLARND( 2, ISEED ) 00176 X( JJ-1 ) = R 00177 JJ = JJ - ( J-1 ) 00178 X( JJ ) = C 00179 JJ = JJ - ( J-2 ) 00180 J = J - 3 00181 END IF 00182 IF( J.GT.1 ) THEN 00183 X( JJ ) = ZLARND( 2, ISEED ) 00184 X( JJ-J ) = ZLARND( 2, ISEED ) 00185 IF( ABS( X( JJ ) ).GT.ABS( X( JJ-J ) ) ) THEN 00186 X( JJ-1 ) = 2.0D0*X( JJ ) 00187 ELSE 00188 X( JJ-1 ) = 2.0D0*X( JJ-J ) 00189 END IF 00190 JJ = JJ - J - ( J-1 ) 00191 J = J - 2 00192 ELSE IF( J.EQ.1 ) THEN 00193 X( JJ ) = ZLARND( 2, ISEED ) 00194 J = J - 1 00195 END IF 00196 * 00197 * UPLO = 'L': Lower triangular storage 00198 * 00199 ELSE 00200 N5 = N / 5 00201 N5 = N5*5 00202 * 00203 JJ = 1 00204 DO 30 J = 1, N5, 5 00205 A = ALPHA3*ZLARND( 5, ISEED ) 00206 B = ZLARND( 5, ISEED ) / ALPHA 00207 C = A - 2.D0*B*EYE 00208 R = C / BETA 00209 X( JJ ) = A 00210 X( JJ+2 ) = B 00211 JJ = JJ + ( N-J+1 ) 00212 X( JJ ) = ZLARND( 2, ISEED ) 00213 X( JJ+1 ) = R 00214 JJ = JJ + ( N-J ) 00215 X( JJ ) = C 00216 JJ = JJ + ( N-J-1 ) 00217 X( JJ ) = ZLARND( 2, ISEED ) 00218 JJ = JJ + ( N-J-2 ) 00219 X( JJ ) = ZLARND( 2, ISEED ) 00220 IF( ABS( X( JJ-( N-J-2 ) ) ).GT.ABS( X( JJ ) ) ) THEN 00221 X( JJ-( N-J-2 )+1 ) = 2.0D0*X( JJ-( N-J-2 ) ) 00222 ELSE 00223 X( JJ-( N-J-2 )+1 ) = 2.0D0*X( JJ ) 00224 END IF 00225 JJ = JJ + ( N-J-3 ) 00226 30 CONTINUE 00227 * 00228 * Clean-up for N not a multiple of 5. 00229 * 00230 J = N5 + 1 00231 IF( J.LT.N-1 ) THEN 00232 A = ALPHA3*ZLARND( 5, ISEED ) 00233 B = ZLARND( 5, ISEED ) / ALPHA 00234 C = A - 2.D0*B*EYE 00235 R = C / BETA 00236 X( JJ ) = A 00237 X( JJ+2 ) = B 00238 JJ = JJ + ( N-J+1 ) 00239 X( JJ ) = ZLARND( 2, ISEED ) 00240 X( JJ+1 ) = R 00241 JJ = JJ + ( N-J ) 00242 X( JJ ) = C 00243 JJ = JJ + ( N-J-1 ) 00244 J = J + 3 00245 END IF 00246 IF( J.LT.N ) THEN 00247 X( JJ ) = ZLARND( 2, ISEED ) 00248 X( JJ+( N-J+1 ) ) = ZLARND( 2, ISEED ) 00249 IF( ABS( X( JJ ) ).GT.ABS( X( JJ+( N-J+1 ) ) ) ) THEN 00250 X( JJ+1 ) = 2.0D0*X( JJ ) 00251 ELSE 00252 X( JJ+1 ) = 2.0D0*X( JJ+( N-J+1 ) ) 00253 END IF 00254 JJ = JJ + ( N-J+1 ) + ( N-J ) 00255 J = J + 2 00256 ELSE IF( J.EQ.N ) THEN 00257 X( JJ ) = ZLARND( 2, ISEED ) 00258 JJ = JJ + ( N-J+1 ) 00259 J = J + 1 00260 END IF 00261 END IF 00262 * 00263 RETURN 00264 * 00265 * End of ZLATSP 00266 * 00267 END