hfti.f 9.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288
  1. *DECK HFTI
  2. SUBROUTINE HFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H,
  3. + G, IP)
  4. C***BEGIN PROLOGUE HFTI
  5. C***PURPOSE Solve a linear least squares problems by performing a QR
  6. C factorization of the matrix using Householder
  7. C transformations.
  8. C***LIBRARY SLATEC
  9. C***CATEGORY D9
  10. C***TYPE SINGLE PRECISION (HFTI-S, DHFTI-D)
  11. C***KEYWORDS CURVE FITTING, LINEAR LEAST SQUARES, QR FACTORIZATION
  12. C***AUTHOR Lawson, C. L., (JPL)
  13. C Hanson, R. J., (SNLA)
  14. C***DESCRIPTION
  15. C
  16. C DIMENSION A(MDA,N),(B(MDB,NB) or B(M)),RNORM(NB),H(N),G(N),IP(N)
  17. C
  18. C This subroutine solves a linear least squares problem or a set of
  19. C linear least squares problems having the same matrix but different
  20. C right-side vectors. The problem data consists of an M by N matrix
  21. C A, an M by NB matrix B, and an absolute tolerance parameter TAU
  22. C whose usage is described below. The NB column vectors of B
  23. C represent right-side vectors for NB distinct linear least squares
  24. C problems.
  25. C
  26. C This set of problems can also be written as the matrix least
  27. C squares problem
  28. C
  29. C AX = B,
  30. C
  31. C where X is the N by NB solution matrix.
  32. C
  33. C Note that if B is the M by M identity matrix, then X will be the
  34. C pseudo-inverse of A.
  35. C
  36. C This subroutine first transforms the augmented matrix (A B) to a
  37. C matrix (R C) using premultiplying Householder transformations with
  38. C column interchanges. All subdiagonal elements in the matrix R are
  39. C zero and its diagonal elements satisfy
  40. C
  41. C ABS(R(I,I)).GE.ABS(R(I+1,I+1)),
  42. C
  43. C I = 1,...,L-1, where
  44. C
  45. C L = MIN(M,N).
  46. C
  47. C The subroutine will compute an integer, KRANK, equal to the number
  48. C of diagonal terms of R that exceed TAU in magnitude. Then a
  49. C solution of minimum Euclidean length is computed using the first
  50. C KRANK rows of (R C).
  51. C
  52. C To be specific we suggest that the user consider an easily
  53. C computable matrix norm, such as, the maximum of all column sums of
  54. C magnitudes.
  55. C
  56. C Now if the relative uncertainty of B is EPS, (norm of uncertainty/
  57. C norm of B), it is suggested that TAU be set approximately equal to
  58. C EPS*(norm of A).
  59. C
  60. C The user must dimension all arrays appearing in the call list..
  61. C A(MDA,N),(B(MDB,NB) or B(M)),RNORM(NB),H(N),G(N),IP(N). This
  62. C permits the solution of a range of problems in the same array
  63. C space.
  64. C
  65. C The entire set of parameters for HFTI are
  66. C
  67. C INPUT..
  68. C
  69. C A(*,*),MDA,M,N The array A(*,*) initially contains the M by N
  70. C matrix A of the least squares problem AX = B.
  71. C The first dimensioning parameter of the array
  72. C A(*,*) is MDA, which must satisfy MDA.GE.M
  73. C Either M.GE.N or M.LT.N is permitted. There
  74. C is no restriction on the rank of A. The
  75. C condition MDA.LT.M is considered an error.
  76. C
  77. C B(*),MDB,NB If NB = 0 the subroutine will perform the
  78. C orthogonal decomposition but will make no
  79. C references to the array B(*). If NB.GT.0
  80. C the array B(*) must initially contain the M by
  81. C NB matrix B of the least squares problem AX =
  82. C B. If NB.GE.2 the array B(*) must be doubly
  83. C subscripted with first dimensioning parameter
  84. C MDB.GE.MAX(M,N). If NB = 1 the array B(*) may
  85. C be either doubly or singly subscripted. In
  86. C the latter case the value of MDB is arbitrary
  87. C but it should be set to some valid integer
  88. C value such as MDB = M.
  89. C
  90. C The condition of NB.GT.1.AND.MDB.LT. MAX(M,N)
  91. C is considered an error.
  92. C
  93. C TAU Absolute tolerance parameter provided by user
  94. C for pseudorank determination.
  95. C
  96. C H(*),G(*),IP(*) Arrays of working space used by HFTI.
  97. C
  98. C OUTPUT..
  99. C
  100. C A(*,*) The contents of the array A(*,*) will be
  101. C modified by the subroutine. These contents
  102. C are not generally required by the user.
  103. C
  104. C B(*) On return the array B(*) will contain the N by
  105. C NB solution matrix X.
  106. C
  107. C KRANK Set by the subroutine to indicate the
  108. C pseudorank of A.
  109. C
  110. C RNORM(*) On return, RNORM(J) will contain the Euclidean
  111. C norm of the residual vector for the problem
  112. C defined by the J-th column vector of the array
  113. C B(*,*) for J = 1,...,NB.
  114. C
  115. C H(*),G(*) On return these arrays respectively contain
  116. C elements of the pre- and post-multiplying
  117. C Householder transformations used to compute
  118. C the minimum Euclidean length solution.
  119. C
  120. C IP(*) Array in which the subroutine records indices
  121. C describing the permutation of column vectors.
  122. C The contents of arrays H(*),G(*) and IP(*)
  123. C are not generally required by the user.
  124. C
  125. C***REFERENCES C. L. Lawson and R. J. Hanson, Solving Least Squares
  126. C Problems, Prentice-Hall, Inc., 1974, Chapter 14.
  127. C***ROUTINES CALLED H12, R1MACH, XERMSG
  128. C***REVISION HISTORY (YYMMDD)
  129. C 790101 DATE WRITTEN
  130. C 890531 Changed all specific intrinsics to generic. (WRB)
  131. C 891006 Cosmetic changes to prologue. (WRB)
  132. C 891006 REVISION DATE from Version 3.2
  133. C 891214 Prologue converted to Version 4.0 format. (BAB)
  134. C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
  135. C 901005 Replace usage of DIFF with usage of R1MACH. (RWC)
  136. C 920501 Reformatted the REFERENCES section. (WRB)
  137. C***END PROLOGUE HFTI
  138. DIMENSION A(MDA,*),B(MDB,*),H(*),G(*),RNORM(*)
  139. INTEGER IP(*)
  140. DOUBLE PRECISION SM,DZERO
  141. SAVE RELEPS
  142. DATA RELEPS /0.E0/
  143. C***FIRST EXECUTABLE STATEMENT HFTI
  144. IF (RELEPS.EQ.0) RELEPS = R1MACH(4)
  145. SZERO=0.
  146. DZERO=0.D0
  147. FACTOR=0.001
  148. C
  149. K=0
  150. LDIAG=MIN(M,N)
  151. IF (LDIAG.LE.0) GO TO 270
  152. IF (.NOT.MDA.LT.M) GO TO 5
  153. NERR=1
  154. IOPT=2
  155. CALL XERMSG ('SLATEC', 'HFTI', 'MDA.LT.M, PROBABLE ERROR.',
  156. + NERR, IOPT)
  157. RETURN
  158. 5 CONTINUE
  159. C
  160. IF (.NOT.(NB.GT.1.AND.MAX(M,N).GT.MDB)) GO TO 6
  161. NERR=2
  162. IOPT=2
  163. CALL XERMSG ('SLATEC', 'HFTI',
  164. + 'MDB.LT.MAX(M,N).AND.NB.GT.1. PROBABLE ERROR.', NERR, IOPT)
  165. RETURN
  166. 6 CONTINUE
  167. C
  168. DO 80 J=1,LDIAG
  169. IF (J.EQ.1) GO TO 20
  170. C
  171. C UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX
  172. C ..
  173. LMAX=J
  174. DO 10 L=J,N
  175. H(L)=H(L)-A(J-1,L)**2
  176. IF (H(L).GT.H(LMAX)) LMAX=L
  177. 10 CONTINUE
  178. IF (FACTOR*H(LMAX) .GT. HMAX*RELEPS) GO TO 50
  179. C
  180. C COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX
  181. C ..
  182. 20 LMAX=J
  183. DO 40 L=J,N
  184. H(L)=0.
  185. DO 30 I=J,M
  186. 30 H(L)=H(L)+A(I,L)**2
  187. IF (H(L).GT.H(LMAX)) LMAX=L
  188. 40 CONTINUE
  189. HMAX=H(LMAX)
  190. C ..
  191. C LMAX HAS BEEN DETERMINED
  192. C
  193. C DO COLUMN INTERCHANGES IF NEEDED.
  194. C ..
  195. 50 CONTINUE
  196. IP(J)=LMAX
  197. IF (IP(J).EQ.J) GO TO 70
  198. DO 60 I=1,M
  199. TMP=A(I,J)
  200. A(I,J)=A(I,LMAX)
  201. 60 A(I,LMAX)=TMP
  202. H(LMAX)=H(J)
  203. C
  204. C COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B.
  205. C ..
  206. 70 CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J)
  207. 80 CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB)
  208. C
  209. C DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU.
  210. C ..
  211. DO 90 J=1,LDIAG
  212. IF (ABS(A(J,J)).LE.TAU) GO TO 100
  213. 90 CONTINUE
  214. K=LDIAG
  215. GO TO 110
  216. 100 K=J-1
  217. 110 KP1=K+1
  218. C
  219. C COMPUTE THE NORMS OF THE RESIDUAL VECTORS.
  220. C
  221. IF (NB.LE.0) GO TO 140
  222. DO 130 JB=1,NB
  223. TMP=SZERO
  224. IF (KP1.GT.M) GO TO 130
  225. DO 120 I=KP1,M
  226. 120 TMP=TMP+B(I,JB)**2
  227. 130 RNORM(JB)=SQRT(TMP)
  228. 140 CONTINUE
  229. C SPECIAL FOR PSEUDORANK = 0
  230. IF (K.GT.0) GO TO 160
  231. IF (NB.LE.0) GO TO 270
  232. DO 150 JB=1,NB
  233. DO 150 I=1,N
  234. 150 B(I,JB)=SZERO
  235. GO TO 270
  236. C
  237. C IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER
  238. C DECOMPOSITION OF FIRST K ROWS.
  239. C ..
  240. 160 IF (K.EQ.N) GO TO 180
  241. DO 170 II=1,K
  242. I=KP1-II
  243. 170 CALL H12 (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1)
  244. 180 CONTINUE
  245. C
  246. C
  247. IF (NB.LE.0) GO TO 270
  248. DO 260 JB=1,NB
  249. C
  250. C SOLVE THE K BY K TRIANGULAR SYSTEM.
  251. C ..
  252. DO 210 L=1,K
  253. SM=DZERO
  254. I=KP1-L
  255. IF (I.EQ.K) GO TO 200
  256. IP1=I+1
  257. DO 190 J=IP1,K
  258. 190 SM=SM+A(I,J)*DBLE(B(J,JB))
  259. 200 SM1=SM
  260. 210 B(I,JB)=(B(I,JB)-SM1)/A(I,I)
  261. C
  262. C COMPLETE COMPUTATION OF SOLUTION VECTOR.
  263. C ..
  264. IF (K.EQ.N) GO TO 240
  265. DO 220 J=KP1,N
  266. 220 B(J,JB)=SZERO
  267. DO 230 I=1,K
  268. 230 CALL H12 (2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1)
  269. C
  270. C RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE
  271. C COLUMN INTERCHANGES.
  272. C ..
  273. 240 DO 250 JJ=1,LDIAG
  274. J=LDIAG+1-JJ
  275. IF (IP(J).EQ.J) GO TO 250
  276. L=IP(J)
  277. TMP=B(L,JB)
  278. B(L,JB)=B(J,JB)
  279. B(J,JB)=TMP
  280. 250 CONTINUE
  281. 260 CONTINUE
  282. C ..
  283. C THE SOLUTION VECTORS, X, ARE NOW
  284. C IN THE FIRST N ROWS OF THE ARRAY B(,).
  285. C
  286. 270 KRANK=K
  287. RETURN
  288. END