00001 C Work performed under the auspices of the U.S. Department of Energy 00002 C by Lawrence Livermore National Laboratory under contract number 00003 C W-7405-Eng-48. 00004 C 00005 SUBROUTINE DHEQR (A, LDA, N, Q, INFO, IJOB) 00006 C 00007 C***BEGIN PROLOGUE DHEQR 00008 C***DATE WRITTEN 890101 (YYMMDD) 00009 C***REVISION DATE 900926 (YYMMDD) 00010 C 00011 C----------------------------------------------------------------------- 00012 C***DESCRIPTION 00013 C 00014 C This routine performs a QR decomposition of an upper 00015 C Hessenberg matrix A. There are two options available: 00016 C 00017 C (1) performing a fresh decomposition 00018 C (2) updating the QR factors by adding a row and A 00019 C column to the matrix A. 00020 C 00021 C DHEQR decomposes an upper Hessenberg matrix by using Givens 00022 C rotations. 00023 C 00024 C On entry 00025 C 00026 C A DOUBLE PRECISION(LDA, N) 00027 C The matrix to be decomposed. 00028 C 00029 C LDA INTEGER 00030 C The leading dimension of the array A. 00031 C 00032 C N INTEGER 00033 C A is an (N+1) by N Hessenberg matrix. 00034 C 00035 C IJOB INTEGER 00036 C = 1 Means that a fresh decomposition of the 00037 C matrix A is desired. 00038 C .GE. 2 Means that the current decomposition of A 00039 C will be updated by the addition of a row 00040 C and a column. 00041 C On return 00042 C 00043 C A The upper triangular matrix R. 00044 C The factorization can be written Q*A = R, where 00045 C Q is a product of Givens rotations and R is upper 00046 C triangular. 00047 C 00048 C Q DOUBLE PRECISION(2*N) 00049 C The factors C and S of each Givens rotation used 00050 C in decomposing A. 00051 C 00052 C INFO INTEGER 00053 C = 0 normal value. 00054 C = K If A(K,K) .EQ. 0.0. This is not an error 00055 C condition for this subroutine, but it does 00056 C indicate that DHELS will divide by zero 00057 C if called. 00058 C 00059 C Modification of LINPACK. 00060 C Peter Brown, Lawrence Livermore Natl. Lab. 00061 C 00062 C----------------------------------------------------------------------- 00063 C***ROUTINES CALLED (NONE) 00064 C 00065 C***END PROLOGUE DHEQR 00066 C 00067 INTEGER LDA, N, INFO, IJOB 00068 DOUBLE PRECISION A(LDA,*), Q(*) 00069 INTEGER I, IQ, J, K, KM1, KP1, NM1 00070 DOUBLE PRECISION C, S, T, T1, T2 00071 C 00072 IF (IJOB .GT. 1) GO TO 70 00073 C----------------------------------------------------------------------- 00074 C A new factorization is desired. 00075 C----------------------------------------------------------------------- 00076 C 00077 C QR decomposition without pivoting. 00078 C 00079 INFO = 0 00080 DO 60 K = 1, N 00081 KM1 = K - 1 00082 KP1 = K + 1 00083 C 00084 C Compute Kth column of R. 00085 C First, multiply the Kth column of A by the previous 00086 C K-1 Givens rotations. 00087 C 00088 IF (KM1 .LT. 1) GO TO 20 00089 DO 10 J = 1, KM1 00090 I = 2*(J-1) + 1 00091 T1 = A(J,K) 00092 T2 = A(J+1,K) 00093 C = Q(I) 00094 S = Q(I+1) 00095 A(J,K) = C*T1 - S*T2 00096 A(J+1,K) = S*T1 + C*T2 00097 10 CONTINUE 00098 C 00099 C Compute Givens components C and S. 00100 C 00101 20 CONTINUE 00102 IQ = 2*KM1 + 1 00103 T1 = A(K,K) 00104 T2 = A(KP1,K) 00105 IF (T2 .NE. 0.0D0) GO TO 30 00106 C = 1.0D0 00107 S = 0.0D0 00108 GO TO 50 00109 30 CONTINUE 00110 IF (ABS(T2) .LT. ABS(T1)) GO TO 40 00111 T = T1/T2 00112 S = -1.0D0/SQRT(1.0D0+T*T) 00113 C = -S*T 00114 GO TO 50 00115 40 CONTINUE 00116 T = T2/T1 00117 C = 1.0D0/SQRT(1.0D0+T*T) 00118 S = -C*T 00119 50 CONTINUE 00120 Q(IQ) = C 00121 Q(IQ+1) = S 00122 A(K,K) = C*T1 - S*T2 00123 IF (A(K,K) .EQ. 0.0D0) INFO = K 00124 60 CONTINUE 00125 RETURN 00126 C----------------------------------------------------------------------- 00127 C The old factorization of A will be updated. A row and a column 00128 C has been added to the matrix A. 00129 C N by N-1 is now the old size of the matrix. 00130 C----------------------------------------------------------------------- 00131 70 CONTINUE 00132 NM1 = N - 1 00133 C----------------------------------------------------------------------- 00134 C Multiply the new column by the N previous Givens rotations. 00135 C----------------------------------------------------------------------- 00136 DO 100 K = 1,NM1 00137 I = 2*(K-1) + 1 00138 T1 = A(K,N) 00139 T2 = A(K+1,N) 00140 C = Q(I) 00141 S = Q(I+1) 00142 A(K,N) = C*T1 - S*T2 00143 A(K+1,N) = S*T1 + C*T2 00144 100 CONTINUE 00145 C----------------------------------------------------------------------- 00146 C Complete update of decomposition by forming last Givens rotation, 00147 C and multiplying it times the column vector (A(N,N),A(NP1,N)). 00148 C----------------------------------------------------------------------- 00149 INFO = 0 00150 T1 = A(N,N) 00151 T2 = A(N+1,N) 00152 IF (T2 .NE. 0.0D0) GO TO 110 00153 C = 1.0D0 00154 S = 0.0D0 00155 GO TO 130 00156 110 CONTINUE 00157 IF (ABS(T2) .LT. ABS(T1)) GO TO 120 00158 T = T1/T2 00159 S = -1.0D0/SQRT(1.0D0+T*T) 00160 C = -S*T 00161 GO TO 130 00162 120 CONTINUE 00163 T = T2/T1 00164 C = 1.0D0/SQRT(1.0D0+T*T) 00165 S = -C*T 00166 130 CONTINUE 00167 IQ = 2*N - 1 00168 Q(IQ) = C 00169 Q(IQ+1) = S 00170 A(N,N) = C*T1 - S*T2 00171 IF (A(N,N) .EQ. 0.0D0) INFO = N 00172 RETURN 00173 C 00174 C------END OF SUBROUTINE DHEQR------------------------------------------ 00175 END
1.7.1
