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 DHELS (A, LDA, N, Q, B) 00006 C 00007 C***BEGIN PROLOGUE DHELS 00008 C***DATE WRITTEN 890101 (YYMMDD) 00009 C***REVISION DATE 900926 (YYMMDD) 00010 C 00011 C 00012 C----------------------------------------------------------------------- 00013 C***DESCRIPTION 00014 C 00015 C This is similar to the LINPACK routine DGESL except that 00016 C A is an upper Hessenberg matrix. 00017 C 00018 C DHELS solves the least squares problem 00019 C 00020 C MIN (B-A*X,B-A*X) 00021 C 00022 C using the factors computed by DHEQR. 00023 C 00024 C On entry 00025 C 00026 C A DOUBLE PRECISION (LDA, N) 00027 C The output from DHEQR which contains the upper 00028 C triangular factor R in the QR decomposition of A. 00029 C 00030 C LDA INTEGER 00031 C The leading dimension of the array A . 00032 C 00033 C N INTEGER 00034 C A is originally an (N+1) by N matrix. 00035 C 00036 C Q DOUBLE PRECISION(2*N) 00037 C The coefficients of the N givens rotations 00038 C used in the QR factorization of A. 00039 C 00040 C B DOUBLE PRECISION(N+1) 00041 C The right hand side vector. 00042 C 00043 C 00044 C On return 00045 C 00046 C B The solution vector X. 00047 C 00048 C 00049 C Modification of LINPACK. 00050 C Peter Brown, Lawrence Livermore Natl. Lab. 00051 C 00052 C----------------------------------------------------------------------- 00053 C***ROUTINES CALLED 00054 C DAXPY 00055 C 00056 C***END PROLOGUE DHELS 00057 C 00058 INTEGER LDA, N 00059 DOUBLE PRECISION A(LDA,*), B(*), Q(*) 00060 INTEGER IQ, K, KB, KP1 00061 DOUBLE PRECISION C, S, T, T1, T2 00062 C 00063 C Minimize (B-A*X,B-A*X). 00064 C First form Q*B. 00065 C 00066 DO 20 K = 1, N 00067 KP1 = K + 1 00068 IQ = 2*(K-1) + 1 00069 C = Q(IQ) 00070 S = Q(IQ+1) 00071 T1 = B(K) 00072 T2 = B(KP1) 00073 B(K) = C*T1 - S*T2 00074 B(KP1) = S*T1 + C*T2 00075 20 CONTINUE 00076 C 00077 C Now solve R*X = Q*B. 00078 C 00079 DO 40 KB = 1, N 00080 K = N + 1 - KB 00081 B(K) = B(K)/A(K,K) 00082 T = -B(K) 00083 CALL DAXPY (K-1, T, A(1,K), 1, B(1), 1) 00084 40 CONTINUE 00085 RETURN 00086 C 00087 C------END OF SUBROUTINE DHELS------------------------------------------ 00088 END
1.7.1
