GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dhels.f
Go to the documentation of this file.
1 C Work performed under the auspices of the U.S. Department of Energy
2 C by Lawrence Livermore National Laboratory under contract number
3 C W-7405-Eng-48.
4 C
5  SUBROUTINE dhels (A, LDA, N, Q, B)
6 C
7 C***BEGIN PROLOGUE DHELS
8 C***DATE WRITTEN 890101 (YYMMDD)
9 C***REVISION DATE 900926 (YYMMDD)
10 C
11 C
12 C-----------------------------------------------------------------------
13 C***DESCRIPTION
14 C
15 C This is similar to the LINPACK routine DGESL except that
16 C A is an upper Hessenberg matrix.
17 C
18 C DHELS solves the least squares problem
19 C
20 C MIN (B-A*X,B-A*X)
21 C
22 C using the factors computed by DHEQR.
23 C
24 C On entry
25 C
26 C A DOUBLE PRECISION (LDA, N)
27 C The output from DHEQR which contains the upper
28 C triangular factor R in the QR decomposition of A.
29 C
30 C LDA INTEGER
31 C The leading dimension of the array A .
32 C
33 C N INTEGER
34 C A is originally an (N+1) by N matrix.
35 C
36 C Q DOUBLE PRECISION(2*N)
37 C The coefficients of the N givens rotations
38 C used in the QR factorization of A.
39 C
40 C B DOUBLE PRECISION(N+1)
41 C The right hand side vector.
42 C
43 C
44 C On return
45 C
46 C B The solution vector X.
47 C
48 C
49 C Modification of LINPACK.
50 C Peter Brown, Lawrence Livermore Natl. Lab.
51 C
52 C-----------------------------------------------------------------------
53 C***ROUTINES CALLED
54 C DAXPY
55 C
56 C***END PROLOGUE DHELS
57 C
58  INTEGER LDA, N
59  DOUBLE PRECISION A(LDA,*), B(*), Q(*)
60  INTEGER IQ, K, KB, KP1
61  DOUBLE PRECISION C, S, T, T1, T2
62 C
63 C Minimize (B-A*X,B-A*X).
64 C First form Q*B.
65 C
66  DO 20 k = 1, n
67  kp1 = k + 1
68  iq = 2*(k-1) + 1
69  c = q(iq)
70  s = q(iq+1)
71  t1 = b(k)
72  t2 = b(kp1)
73  b(k) = c*t1 - s*t2
74  b(kp1) = s*t1 + c*t2
75  20 CONTINUE
76 C
77 C Now solve R*X = Q*B.
78 C
79  DO 40 kb = 1, n
80  k = n + 1 - kb
81  b(k) = b(k)/a(k,k)
82  t = -b(k)
83  CALL daxpy (k-1, t, a(1,k), 1, b(1), 1)
84  40 CONTINUE
85  RETURN
86 C
87 C------END OF SUBROUTINE DHELS------------------------------------------
88  END
subroutine dhels(A, LDA, N, Q, B)
Definition: dhels.f:6