GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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