LAPACK 3.3.1
Linear Algebra PACKage
|
00001 REAL FUNCTION SLANSY( NORM, UPLO, N, A, LDA, WORK ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER NORM, UPLO 00010 INTEGER LDA, N 00011 * .. 00012 * .. Array Arguments .. 00013 REAL A( LDA, * ), WORK( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * SLANSY returns the value of the one norm, or the Frobenius norm, or 00020 * the infinity norm, or the element of largest absolute value of a 00021 * real symmetric matrix A. 00022 * 00023 * Description 00024 * =========== 00025 * 00026 * SLANSY returns the value 00027 * 00028 * SLANSY = ( max(abs(A(i,j))), NORM = 'M' or 'm' 00029 * ( 00030 * ( norm1(A), NORM = '1', 'O' or 'o' 00031 * ( 00032 * ( normI(A), NORM = 'I' or 'i' 00033 * ( 00034 * ( normF(A), NORM = 'F', 'f', 'E' or 'e' 00035 * 00036 * where norm1 denotes the one norm of a matrix (maximum column sum), 00037 * normI denotes the infinity norm of a matrix (maximum row sum) and 00038 * normF denotes the Frobenius norm of a matrix (square root of sum of 00039 * squares). Note that max(abs(A(i,j))) is not a consistent matrix norm. 00040 * 00041 * Arguments 00042 * ========= 00043 * 00044 * NORM (input) CHARACTER*1 00045 * Specifies the value to be returned in SLANSY as described 00046 * above. 00047 * 00048 * UPLO (input) CHARACTER*1 00049 * Specifies whether the upper or lower triangular part of the 00050 * symmetric matrix A is to be referenced. 00051 * = 'U': Upper triangular part of A is referenced 00052 * = 'L': Lower triangular part of A is referenced 00053 * 00054 * N (input) INTEGER 00055 * The order of the matrix A. N >= 0. When N = 0, SLANSY is 00056 * set to zero. 00057 * 00058 * A (input) REAL array, dimension (LDA,N) 00059 * The symmetric matrix A. If UPLO = 'U', the leading n by n 00060 * upper triangular part of A contains the upper triangular part 00061 * of the matrix A, and the strictly lower triangular part of A 00062 * is not referenced. If UPLO = 'L', the leading n by n lower 00063 * triangular part of A contains the lower triangular part of 00064 * the matrix A, and the strictly upper triangular part of A is 00065 * not referenced. 00066 * 00067 * LDA (input) INTEGER 00068 * The leading dimension of the array A. LDA >= max(N,1). 00069 * 00070 * WORK (workspace) REAL array, dimension (MAX(1,LWORK)), 00071 * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, 00072 * WORK is not referenced. 00073 * 00074 * ===================================================================== 00075 * 00076 * .. Parameters .. 00077 REAL ONE, ZERO 00078 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00079 * .. 00080 * .. Local Scalars .. 00081 INTEGER I, J 00082 REAL ABSA, SCALE, SUM, VALUE 00083 * .. 00084 * .. External Subroutines .. 00085 EXTERNAL SLASSQ 00086 * .. 00087 * .. External Functions .. 00088 LOGICAL LSAME 00089 EXTERNAL LSAME 00090 * .. 00091 * .. Intrinsic Functions .. 00092 INTRINSIC ABS, MAX, SQRT 00093 * .. 00094 * .. Executable Statements .. 00095 * 00096 IF( N.EQ.0 ) THEN 00097 VALUE = ZERO 00098 ELSE IF( LSAME( NORM, 'M' ) ) THEN 00099 * 00100 * Find max(abs(A(i,j))). 00101 * 00102 VALUE = ZERO 00103 IF( LSAME( UPLO, 'U' ) ) THEN 00104 DO 20 J = 1, N 00105 DO 10 I = 1, J 00106 VALUE = MAX( VALUE, ABS( A( I, J ) ) ) 00107 10 CONTINUE 00108 20 CONTINUE 00109 ELSE 00110 DO 40 J = 1, N 00111 DO 30 I = J, N 00112 VALUE = MAX( VALUE, ABS( A( I, J ) ) ) 00113 30 CONTINUE 00114 40 CONTINUE 00115 END IF 00116 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. 00117 $ ( NORM.EQ.'1' ) ) THEN 00118 * 00119 * Find normI(A) ( = norm1(A), since A is symmetric). 00120 * 00121 VALUE = ZERO 00122 IF( LSAME( UPLO, 'U' ) ) THEN 00123 DO 60 J = 1, N 00124 SUM = ZERO 00125 DO 50 I = 1, J - 1 00126 ABSA = ABS( A( I, J ) ) 00127 SUM = SUM + ABSA 00128 WORK( I ) = WORK( I ) + ABSA 00129 50 CONTINUE 00130 WORK( J ) = SUM + ABS( A( J, J ) ) 00131 60 CONTINUE 00132 DO 70 I = 1, N 00133 VALUE = MAX( VALUE, WORK( I ) ) 00134 70 CONTINUE 00135 ELSE 00136 DO 80 I = 1, N 00137 WORK( I ) = ZERO 00138 80 CONTINUE 00139 DO 100 J = 1, N 00140 SUM = WORK( J ) + ABS( A( J, J ) ) 00141 DO 90 I = J + 1, N 00142 ABSA = ABS( A( I, J ) ) 00143 SUM = SUM + ABSA 00144 WORK( I ) = WORK( I ) + ABSA 00145 90 CONTINUE 00146 VALUE = MAX( VALUE, SUM ) 00147 100 CONTINUE 00148 END IF 00149 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 00150 * 00151 * Find normF(A). 00152 * 00153 SCALE = ZERO 00154 SUM = ONE 00155 IF( LSAME( UPLO, 'U' ) ) THEN 00156 DO 110 J = 2, N 00157 CALL SLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) 00158 110 CONTINUE 00159 ELSE 00160 DO 120 J = 1, N - 1 00161 CALL SLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) 00162 120 CONTINUE 00163 END IF 00164 SUM = 2*SUM 00165 CALL SLASSQ( N, A, LDA+1, SCALE, SUM ) 00166 VALUE = SCALE*SQRT( SUM ) 00167 END IF 00168 * 00169 SLANSY = VALUE 00170 RETURN 00171 * 00172 * End of SLANSY 00173 * 00174 END