dlsei.f 28 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735
  1. *DECK DLSEI
  2. SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME,
  3. + RNORML, MODE, WS, IP)
  4. C***BEGIN PROLOGUE DLSEI
  5. C***PURPOSE Solve a linearly constrained least squares problem with
  6. C equality and inequality constraints, and optionally compute
  7. C a covariance matrix.
  8. C***LIBRARY SLATEC
  9. C***CATEGORY K1A2A, D9
  10. C***TYPE DOUBLE PRECISION (LSEI-S, DLSEI-D)
  11. C***KEYWORDS CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING,
  12. C EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS,
  13. C QUADRATIC PROGRAMMING
  14. C***AUTHOR Hanson, R. J., (SNLA)
  15. C Haskell, K. H., (SNLA)
  16. C***DESCRIPTION
  17. C
  18. C Abstract
  19. C
  20. C This subprogram solves a linearly constrained least squares
  21. C problem with both equality and inequality constraints, and, if the
  22. C user requests, obtains a covariance matrix of the solution
  23. C parameters.
  24. C
  25. C Suppose there are given matrices E, A and G of respective
  26. C dimensions ME by N, MA by N and MG by N, and vectors F, B and H of
  27. C respective lengths ME, MA and MG. This subroutine solves the
  28. C linearly constrained least squares problem
  29. C
  30. C EX = F, (E ME by N) (equations to be exactly
  31. C satisfied)
  32. C AX = B, (A MA by N) (equations to be
  33. C approximately satisfied,
  34. C least squares sense)
  35. C GX .GE. H,(G MG by N) (inequality constraints)
  36. C
  37. C The inequalities GX .GE. H mean that every component of the
  38. C product GX must be .GE. the corresponding component of H.
  39. C
  40. C In case the equality constraints cannot be satisfied, a
  41. C generalized inverse solution residual vector length is obtained
  42. C for F-EX. This is the minimal length possible for F-EX.
  43. C
  44. C Any values ME .GE. 0, MA .GE. 0, or MG .GE. 0 are permitted. The
  45. C rank of the matrix E is estimated during the computation. We call
  46. C this value KRANKE. It is an output parameter in IP(1) defined
  47. C below. Using a generalized inverse solution of EX=F, a reduced
  48. C least squares problem with inequality constraints is obtained.
  49. C The tolerances used in these tests for determining the rank
  50. C of E and the rank of the reduced least squares problem are
  51. C given in Sandia Tech. Rept. SAND-78-1290. They can be
  52. C modified by the user if new values are provided in
  53. C the option list of the array PRGOPT(*).
  54. C
  55. C The user must dimension all arrays appearing in the call list..
  56. C W(MDW,N+1),PRGOPT(*),X(N),WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2)
  57. C where K=MAX(MA+MG,N). This allows for a solution of a range of
  58. C problems in the given working space. The dimension of WS(*)
  59. C given is a necessary overestimate. Once a particular problem
  60. C has been run, the output parameter IP(3) gives the actual
  61. C dimension required for that problem.
  62. C
  63. C The parameters for DLSEI( ) are
  64. C
  65. C Input.. All TYPE REAL variables are DOUBLE PRECISION
  66. C
  67. C W(*,*),MDW, The array W(*,*) is doubly subscripted with
  68. C ME,MA,MG,N first dimensioning parameter equal to MDW.
  69. C For this discussion let us call M = ME+MA+MG. Then
  70. C MDW must satisfy MDW .GE. M. The condition
  71. C MDW .LT. M is an error.
  72. C
  73. C The array W(*,*) contains the matrices and vectors
  74. C
  75. C (E F)
  76. C (A B)
  77. C (G H)
  78. C
  79. C in rows and columns 1,...,M and 1,...,N+1
  80. C respectively.
  81. C
  82. C The integers ME, MA, and MG are the
  83. C respective matrix row dimensions
  84. C of E, A and G. Each matrix has N columns.
  85. C
  86. C PRGOPT(*) This real-valued array is the option vector.
  87. C If the user is satisfied with the nominal
  88. C subprogram features set
  89. C
  90. C PRGOPT(1)=1 (or PRGOPT(1)=1.0)
  91. C
  92. C Otherwise PRGOPT(*) is a linked list consisting of
  93. C groups of data of the following form
  94. C
  95. C LINK
  96. C KEY
  97. C DATA SET
  98. C
  99. C The parameters LINK and KEY are each one word.
  100. C The DATA SET can be comprised of several words.
  101. C The number of items depends on the value of KEY.
  102. C The value of LINK points to the first
  103. C entry of the next group of data within
  104. C PRGOPT(*). The exception is when there are
  105. C no more options to change. In that
  106. C case, LINK=1 and the values KEY and DATA SET
  107. C are not referenced. The general layout of
  108. C PRGOPT(*) is as follows.
  109. C
  110. C ...PRGOPT(1) = LINK1 (link to first entry of next group)
  111. C . PRGOPT(2) = KEY1 (key to the option change)
  112. C . PRGOPT(3) = data value (data value for this change)
  113. C . .
  114. C . .
  115. C . .
  116. C ...PRGOPT(LINK1) = LINK2 (link to the first entry of
  117. C . next group)
  118. C . PRGOPT(LINK1+1) = KEY2 (key to the option change)
  119. C . PRGOPT(LINK1+2) = data value
  120. C ... .
  121. C . .
  122. C . .
  123. C ...PRGOPT(LINK) = 1 (no more options to change)
  124. C
  125. C Values of LINK that are nonpositive are errors.
  126. C A value of LINK .GT. NLINK=100000 is also an error.
  127. C This helps prevent using invalid but positive
  128. C values of LINK that will probably extend
  129. C beyond the program limits of PRGOPT(*).
  130. C Unrecognized values of KEY are ignored. The
  131. C order of the options is arbitrary and any number
  132. C of options can be changed with the following
  133. C restriction. To prevent cycling in the
  134. C processing of the option array, a count of the
  135. C number of options changed is maintained.
  136. C Whenever this count exceeds NOPT=1000, an error
  137. C message is printed and the subprogram returns.
  138. C
  139. C Options..
  140. C
  141. C KEY=1
  142. C Compute in W(*,*) the N by N
  143. C covariance matrix of the solution variables
  144. C as an output parameter. Nominally the
  145. C covariance matrix will not be computed.
  146. C (This requires no user input.)
  147. C The data set for this option is a single value.
  148. C It must be nonzero when the covariance matrix
  149. C is desired. If it is zero, the covariance
  150. C matrix is not computed. When the covariance matrix
  151. C is computed, the first dimensioning parameter
  152. C of the array W(*,*) must satisfy MDW .GE. MAX(M,N).
  153. C
  154. C KEY=10
  155. C Suppress scaling of the inverse of the
  156. C normal matrix by the scale factor RNORM**2/
  157. C MAX(1, no. of degrees of freedom). This option
  158. C only applies when the option for computing the
  159. C covariance matrix (KEY=1) is used. With KEY=1 and
  160. C KEY=10 used as options the unscaled inverse of the
  161. C normal matrix is returned in W(*,*).
  162. C The data set for this option is a single value.
  163. C When it is nonzero no scaling is done. When it is
  164. C zero scaling is done. The nominal case is to do
  165. C scaling so if option (KEY=1) is used alone, the
  166. C matrix will be scaled on output.
  167. C
  168. C KEY=2
  169. C Scale the nonzero columns of the
  170. C entire data matrix.
  171. C (E)
  172. C (A)
  173. C (G)
  174. C
  175. C to have length one. The data set for this
  176. C option is a single value. It must be
  177. C nonzero if unit length column scaling
  178. C is desired.
  179. C
  180. C KEY=3
  181. C Scale columns of the entire data matrix
  182. C (E)
  183. C (A)
  184. C (G)
  185. C
  186. C with a user-provided diagonal matrix.
  187. C The data set for this option consists
  188. C of the N diagonal scaling factors, one for
  189. C each matrix column.
  190. C
  191. C KEY=4
  192. C Change the rank determination tolerance for
  193. C the equality constraint equations from
  194. C the nominal value of SQRT(DRELPR). This quantity can
  195. C be no smaller than DRELPR, the arithmetic-
  196. C storage precision. The quantity DRELPR is the
  197. C largest positive number such that T=1.+DRELPR
  198. C satisfies T .EQ. 1. The quantity used
  199. C here is internally restricted to be at
  200. C least DRELPR. The data set for this option
  201. C is the new tolerance.
  202. C
  203. C KEY=5
  204. C Change the rank determination tolerance for
  205. C the reduced least squares equations from
  206. C the nominal value of SQRT(DRELPR). This quantity can
  207. C be no smaller than DRELPR, the arithmetic-
  208. C storage precision. The quantity used
  209. C here is internally restricted to be at
  210. C least DRELPR. The data set for this option
  211. C is the new tolerance.
  212. C
  213. C For example, suppose we want to change
  214. C the tolerance for the reduced least squares
  215. C problem, compute the covariance matrix of
  216. C the solution parameters, and provide
  217. C column scaling for the data matrix. For
  218. C these options the dimension of PRGOPT(*)
  219. C must be at least N+9. The Fortran statements
  220. C defining these options would be as follows:
  221. C
  222. C PRGOPT(1)=4 (link to entry 4 in PRGOPT(*))
  223. C PRGOPT(2)=1 (covariance matrix key)
  224. C PRGOPT(3)=1 (covariance matrix wanted)
  225. C
  226. C PRGOPT(4)=7 (link to entry 7 in PRGOPT(*))
  227. C PRGOPT(5)=5 (least squares equas. tolerance key)
  228. C PRGOPT(6)=... (new value of the tolerance)
  229. C
  230. C PRGOPT(7)=N+9 (link to entry N+9 in PRGOPT(*))
  231. C PRGOPT(8)=3 (user-provided column scaling key)
  232. C
  233. C CALL DCOPY (N, D, 1, PRGOPT(9), 1) (Copy the N
  234. C scaling factors from the user array D(*)
  235. C to PRGOPT(9)-PRGOPT(N+8))
  236. C
  237. C PRGOPT(N+9)=1 (no more options to change)
  238. C
  239. C The contents of PRGOPT(*) are not modified
  240. C by the subprogram.
  241. C The options for WNNLS( ) can also be included
  242. C in this array. The values of KEY recognized
  243. C by WNNLS( ) are 6, 7 and 8. Their functions
  244. C are documented in the usage instructions for
  245. C subroutine WNNLS( ). Normally these options
  246. C do not need to be modified when using DLSEI( ).
  247. C
  248. C IP(1), The amounts of working storage actually
  249. C IP(2) allocated for the working arrays WS(*) and
  250. C IP(*), respectively. These quantities are
  251. C compared with the actual amounts of storage
  252. C needed by DLSEI( ). Insufficient storage
  253. C allocated for either WS(*) or IP(*) is an
  254. C error. This feature was included in DLSEI( )
  255. C because miscalculating the storage formulas
  256. C for WS(*) and IP(*) might very well lead to
  257. C subtle and hard-to-find execution errors.
  258. C
  259. C The length of WS(*) must be at least
  260. C
  261. C LW = 2*(ME+N)+K+(MG+2)*(N+7)
  262. C
  263. C where K = max(MA+MG,N)
  264. C This test will not be made if IP(1).LE.0.
  265. C
  266. C The length of IP(*) must be at least
  267. C
  268. C LIP = MG+2*N+2
  269. C This test will not be made if IP(2).LE.0.
  270. C
  271. C Output.. All TYPE REAL variables are DOUBLE PRECISION
  272. C
  273. C X(*),RNORME, The array X(*) contains the solution parameters
  274. C RNORML if the integer output flag MODE = 0 or 1.
  275. C The definition of MODE is given directly below.
  276. C When MODE = 0 or 1, RNORME and RNORML
  277. C respectively contain the residual vector
  278. C Euclidean lengths of F - EX and B - AX. When
  279. C MODE=1 the equality constraint equations EX=F
  280. C are contradictory, so RNORME .NE. 0. The residual
  281. C vector F-EX has minimal Euclidean length. For
  282. C MODE .GE. 2, none of these parameters is defined.
  283. C
  284. C MODE Integer flag that indicates the subprogram
  285. C status after completion. If MODE .GE. 2, no
  286. C solution has been computed.
  287. C
  288. C MODE =
  289. C
  290. C 0 Both equality and inequality constraints
  291. C are compatible and have been satisfied.
  292. C
  293. C 1 Equality constraints are contradictory.
  294. C A generalized inverse solution of EX=F was used
  295. C to minimize the residual vector length F-EX.
  296. C In this sense, the solution is still meaningful.
  297. C
  298. C 2 Inequality constraints are contradictory.
  299. C
  300. C 3 Both equality and inequality constraints
  301. C are contradictory.
  302. C
  303. C The following interpretation of
  304. C MODE=1,2 or 3 must be made. The
  305. C sets consisting of all solutions
  306. C of the equality constraints EX=F
  307. C and all vectors satisfying GX .GE. H
  308. C have no points in common. (In
  309. C particular this does not say that
  310. C each individual set has no points
  311. C at all, although this could be the
  312. C case.)
  313. C
  314. C 4 Usage error occurred. The value
  315. C of MDW is .LT. ME+MA+MG, MDW is
  316. C .LT. N and a covariance matrix is
  317. C requested, or the option vector
  318. C PRGOPT(*) is not properly defined,
  319. C or the lengths of the working arrays
  320. C WS(*) and IP(*), when specified in
  321. C IP(1) and IP(2) respectively, are not
  322. C long enough.
  323. C
  324. C W(*,*) The array W(*,*) contains the N by N symmetric
  325. C covariance matrix of the solution parameters,
  326. C provided this was requested on input with
  327. C the option vector PRGOPT(*) and the output
  328. C flag is returned with MODE = 0 or 1.
  329. C
  330. C IP(*) The integer working array has three entries
  331. C that provide rank and working array length
  332. C information after completion.
  333. C
  334. C IP(1) = rank of equality constraint
  335. C matrix. Define this quantity
  336. C as KRANKE.
  337. C
  338. C IP(2) = rank of reduced least squares
  339. C problem.
  340. C
  341. C IP(3) = the amount of storage in the
  342. C working array WS(*) that was
  343. C actually used by the subprogram.
  344. C The formula given above for the length
  345. C of WS(*) is a necessary overestimate.
  346. C If exactly the same problem matrices
  347. C are used in subsequent executions,
  348. C the declared dimension of WS(*) can
  349. C be reduced to this output value.
  350. C User Designated
  351. C Working Arrays..
  352. C
  353. C WS(*),IP(*) These are respectively type real
  354. C and type integer working arrays.
  355. C Their required minimal lengths are
  356. C given above.
  357. C
  358. C***REFERENCES K. H. Haskell and R. J. Hanson, An algorithm for
  359. C linear least squares problems with equality and
  360. C nonnegativity constraints, Report SAND77-0552, Sandia
  361. C Laboratories, June 1978.
  362. C K. H. Haskell and R. J. Hanson, Selected algorithms for
  363. C the linearly constrained least squares problem - a
  364. C users guide, Report SAND78-1290, Sandia Laboratories,
  365. C August 1979.
  366. C K. H. Haskell and R. J. Hanson, An algorithm for
  367. C linear least squares problems with equality and
  368. C nonnegativity constraints, Mathematical Programming
  369. C 21 (1981), pp. 98-118.
  370. C R. J. Hanson and K. H. Haskell, Two algorithms for the
  371. C linearly constrained least squares problem, ACM
  372. C Transactions on Mathematical Software, September 1982.
  373. C***ROUTINES CALLED D1MACH, DASUM, DAXPY, DCOPY, DDOT, DH12, DLSI,
  374. C DNRM2, DSCAL, DSWAP, XERMSG
  375. C***REVISION HISTORY (YYMMDD)
  376. C 790701 DATE WRITTEN
  377. C 890531 Changed all specific intrinsics to generic. (WRB)
  378. C 890618 Completely restructured and extensively revised (WRB & RWC)
  379. C 890831 REVISION DATE from Version 3.2
  380. C 891214 Prologue converted to Version 4.0 format. (BAB)
  381. C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
  382. C 900510 Convert XERRWV calls to XERMSG calls. (RWC)
  383. C 900604 DP version created from SP version. (RWC)
  384. C 920501 Reformatted the REFERENCES section. (WRB)
  385. C***END PROLOGUE DLSEI
  386. INTEGER IP(3), MA, MDW, ME, MG, MODE, N
  387. DOUBLE PRECISION PRGOPT(*), RNORME, RNORML, W(MDW,*), WS(*), X(*)
  388. C
  389. EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DDOT, DH12, DLSI, DNRM2,
  390. * DSCAL, DSWAP, XERMSG
  391. DOUBLE PRECISION D1MACH, DASUM, DDOT, DNRM2
  392. C
  393. DOUBLE PRECISION DRELPR, ENORM, FNORM, GAM, RB, RN, RNMAX, SIZE,
  394. * SN, SNMAX, T, TAU, UJ, UP, VJ, XNORM, XNRME
  395. INTEGER I, IMAX, J, JP1, K, KEY, KRANKE, LAST, LCHK, LINK, M,
  396. * MAPKE1, MDEQC, MEND, MEP1, N1, N2, NEXT, NLINK, NOPT, NP1,
  397. * NTIMES
  398. LOGICAL COV, FIRST
  399. CHARACTER*8 XERN1, XERN2, XERN3, XERN4
  400. SAVE FIRST, DRELPR
  401. C
  402. DATA FIRST /.TRUE./
  403. C***FIRST EXECUTABLE STATEMENT DLSEI
  404. C
  405. C Set the nominal tolerance used in the code for the equality
  406. C constraint equations.
  407. C
  408. IF (FIRST) DRELPR = D1MACH(4)
  409. FIRST = .FALSE.
  410. TAU = SQRT(DRELPR)
  411. C
  412. C Check that enough storage was allocated in WS(*) and IP(*).
  413. C
  414. MODE = 4
  415. IF (MIN(N,ME,MA,MG) .LT. 0) THEN
  416. WRITE (XERN1, '(I8)') N
  417. WRITE (XERN2, '(I8)') ME
  418. WRITE (XERN3, '(I8)') MA
  419. WRITE (XERN4, '(I8)') MG
  420. CALL XERMSG ('SLATEC', 'LSEI', 'ALL OF THE VARIABLES N, ME,' //
  421. * ' MA, MG MUST BE .GE. 0$$ENTERED ROUTINE WITH' //
  422. * '$$N = ' // XERN1 //
  423. * '$$ME = ' // XERN2 //
  424. * '$$MA = ' // XERN3 //
  425. * '$$MG = ' // XERN4, 2, 1)
  426. RETURN
  427. ENDIF
  428. C
  429. IF (IP(1).GT.0) THEN
  430. LCHK = 2*(ME+N) + MAX(MA+MG,N) + (MG+2)*(N+7)
  431. IF (IP(1).LT.LCHK) THEN
  432. WRITE (XERN1, '(I8)') LCHK
  433. CALL XERMSG ('SLATEC', 'DLSEI', 'INSUFFICIENT STORAGE ' //
  434. * 'ALLOCATED FOR WS(*), NEED LW = ' // XERN1, 2, 1)
  435. RETURN
  436. ENDIF
  437. ENDIF
  438. C
  439. IF (IP(2).GT.0) THEN
  440. LCHK = MG + 2*N + 2
  441. IF (IP(2).LT.LCHK) THEN
  442. WRITE (XERN1, '(I8)') LCHK
  443. CALL XERMSG ('SLATEC', 'DLSEI', 'INSUFFICIENT STORAGE ' //
  444. * 'ALLOCATED FOR IP(*), NEED LIP = ' // XERN1, 2, 1)
  445. RETURN
  446. ENDIF
  447. ENDIF
  448. C
  449. C Compute number of possible right multiplying Householder
  450. C transformations.
  451. C
  452. M = ME + MA + MG
  453. IF (N.LE.0 .OR. M.LE.0) THEN
  454. MODE = 0
  455. RNORME = 0
  456. RNORML = 0
  457. RETURN
  458. ENDIF
  459. C
  460. IF (MDW.LT.M) THEN
  461. CALL XERMSG ('SLATEC', 'DLSEI', 'MDW.LT.ME+MA+MG IS AN ERROR',
  462. + 2, 1)
  463. RETURN
  464. ENDIF
  465. C
  466. NP1 = N + 1
  467. KRANKE = MIN(ME,N)
  468. N1 = 2*KRANKE + 1
  469. N2 = N1 + N
  470. C
  471. C Set nominal values.
  472. C
  473. C The nominal column scaling used in the code is
  474. C the identity scaling.
  475. C
  476. CALL DCOPY (N, 1.D0, 0, WS(N1), 1)
  477. C
  478. C No covariance matrix is nominally computed.
  479. C
  480. COV = .FALSE.
  481. C
  482. C Process option vector.
  483. C Define bound for number of options to change.
  484. C
  485. NOPT = 1000
  486. NTIMES = 0
  487. C
  488. C Define bound for positive values of LINK.
  489. C
  490. NLINK = 100000
  491. LAST = 1
  492. LINK = PRGOPT(1)
  493. IF (LINK.EQ.0 .OR. LINK.GT.NLINK) THEN
  494. CALL XERMSG ('SLATEC', 'DLSEI',
  495. + 'THE OPTION VECTOR IS UNDEFINED', 2, 1)
  496. RETURN
  497. ENDIF
  498. C
  499. 100 IF (LINK.GT.1) THEN
  500. NTIMES = NTIMES + 1
  501. IF (NTIMES.GT.NOPT) THEN
  502. CALL XERMSG ('SLATEC', 'DLSEI',
  503. + 'THE LINKS IN THE OPTION VECTOR ARE CYCLING.', 2, 1)
  504. RETURN
  505. ENDIF
  506. C
  507. KEY = PRGOPT(LAST+1)
  508. IF (KEY.EQ.1) THEN
  509. COV = PRGOPT(LAST+2) .NE. 0.D0
  510. ELSEIF (KEY.EQ.2 .AND. PRGOPT(LAST+2).NE.0.D0) THEN
  511. DO 110 J = 1,N
  512. T = DNRM2(M,W(1,J),1)
  513. IF (T.NE.0.D0) T = 1.D0/T
  514. WS(J+N1-1) = T
  515. 110 CONTINUE
  516. ELSEIF (KEY.EQ.3) THEN
  517. CALL DCOPY (N, PRGOPT(LAST+2), 1, WS(N1), 1)
  518. ELSEIF (KEY.EQ.4) THEN
  519. TAU = MAX(DRELPR,PRGOPT(LAST+2))
  520. ENDIF
  521. C
  522. NEXT = PRGOPT(LINK)
  523. IF (NEXT.LE.0 .OR. NEXT.GT.NLINK) THEN
  524. CALL XERMSG ('SLATEC', 'DLSEI',
  525. + 'THE OPTION VECTOR IS UNDEFINED', 2, 1)
  526. RETURN
  527. ENDIF
  528. C
  529. LAST = LINK
  530. LINK = NEXT
  531. GO TO 100
  532. ENDIF
  533. C
  534. DO 120 J = 1,N
  535. CALL DSCAL (M, WS(N1+J-1), W(1,J), 1)
  536. 120 CONTINUE
  537. C
  538. IF (COV .AND. MDW.LT.N) THEN
  539. CALL XERMSG ('SLATEC', 'DLSEI',
  540. + 'MDW .LT. N WHEN COV MATRIX NEEDED, IS AN ERROR', 2, 1)
  541. RETURN
  542. ENDIF
  543. C
  544. C Problem definition and option vector OK.
  545. C
  546. MODE = 0
  547. C
  548. C Compute norm of equality constraint matrix and right side.
  549. C
  550. ENORM = 0.D0
  551. DO 130 J = 1,N
  552. ENORM = MAX(ENORM,DASUM(ME,W(1,J),1))
  553. 130 CONTINUE
  554. C
  555. FNORM = DASUM(ME,W(1,NP1),1)
  556. SNMAX = 0.D0
  557. RNMAX = 0.D0
  558. DO 150 I = 1,KRANKE
  559. C
  560. C Compute maximum ratio of vector lengths. Partition is at
  561. C column I.
  562. C
  563. DO 140 K = I,ME
  564. SN = DDOT(N-I+1,W(K,I),MDW,W(K,I),MDW)
  565. RN = DDOT(I-1,W(K,1),MDW,W(K,1),MDW)
  566. IF (RN.EQ.0.D0 .AND. SN.GT.SNMAX) THEN
  567. SNMAX = SN
  568. IMAX = K
  569. ELSEIF (K.EQ.I .OR. SN*RNMAX.GT.RN*SNMAX) THEN
  570. SNMAX = SN
  571. RNMAX = RN
  572. IMAX = K
  573. ENDIF
  574. 140 CONTINUE
  575. C
  576. C Interchange rows if necessary.
  577. C
  578. IF (I.NE.IMAX) CALL DSWAP (NP1, W(I,1), MDW, W(IMAX,1), MDW)
  579. IF (SNMAX.GT.RNMAX*TAU**2) THEN
  580. C
  581. C Eliminate elements I+1,...,N in row I.
  582. C
  583. CALL DH12 (1, I, I+1, N, W(I,1), MDW, WS(I), W(I+1,1), MDW,
  584. + 1, M-I)
  585. ELSE
  586. KRANKE = I - 1
  587. GO TO 160
  588. ENDIF
  589. 150 CONTINUE
  590. C
  591. C Save diagonal terms of lower trapezoidal matrix.
  592. C
  593. 160 CALL DCOPY (KRANKE, W, MDW+1, WS(KRANKE+1), 1)
  594. C
  595. C Use Householder transformation from left to achieve
  596. C KRANKE by KRANKE upper triangular form.
  597. C
  598. IF (KRANKE.LT.ME) THEN
  599. DO 170 K = KRANKE,1,-1
  600. C
  601. C Apply transformation to matrix cols. 1,...,K-1.
  602. C
  603. CALL DH12 (1, K, KRANKE+1, ME, W(1,K), 1, UP, W, 1, MDW,
  604. * K-1)
  605. C
  606. C Apply to rt side vector.
  607. C
  608. CALL DH12 (2, K, KRANKE+1, ME, W(1,K), 1, UP, W(1,NP1), 1,
  609. + 1, 1)
  610. 170 CONTINUE
  611. ENDIF
  612. C
  613. C Solve for variables 1,...,KRANKE in new coordinates.
  614. C
  615. CALL DCOPY (KRANKE, W(1, NP1), 1, X, 1)
  616. DO 180 I = 1,KRANKE
  617. X(I) = (X(I)-DDOT(I-1,W(I,1),MDW,X,1))/W(I,I)
  618. 180 CONTINUE
  619. C
  620. C Compute residuals for reduced problem.
  621. C
  622. MEP1 = ME + 1
  623. RNORML = 0.D0
  624. DO 190 I = MEP1,M
  625. W(I,NP1) = W(I,NP1) - DDOT(KRANKE,W(I,1),MDW,X,1)
  626. SN = DDOT(KRANKE,W(I,1),MDW,W(I,1),MDW)
  627. RN = DDOT(N-KRANKE,W(I,KRANKE+1),MDW,W(I,KRANKE+1),MDW)
  628. IF (RN.LE.SN*TAU**2 .AND. KRANKE.LT.N)
  629. * CALL DCOPY (N-KRANKE, 0.D0, 0, W(I,KRANKE+1), MDW)
  630. 190 CONTINUE
  631. C
  632. C Compute equality constraint equations residual length.
  633. C
  634. RNORME = DNRM2(ME-KRANKE,W(KRANKE+1,NP1),1)
  635. C
  636. C Move reduced problem data upward if KRANKE.LT.ME.
  637. C
  638. IF (KRANKE.LT.ME) THEN
  639. DO 200 J = 1,NP1
  640. CALL DCOPY (M-ME, W(ME+1,J), 1, W(KRANKE+1,J), 1)
  641. 200 CONTINUE
  642. ENDIF
  643. C
  644. C Compute solution of reduced problem.
  645. C
  646. CALL DLSI(W(KRANKE+1, KRANKE+1), MDW, MA, MG, N-KRANKE, PRGOPT,
  647. + X(KRANKE+1), RNORML, MODE, WS(N2), IP(2))
  648. C
  649. C Test for consistency of equality constraints.
  650. C
  651. IF (ME.GT.0) THEN
  652. MDEQC = 0
  653. XNRME = DASUM(KRANKE,W(1,NP1),1)
  654. IF (RNORME.GT.TAU*(ENORM*XNRME+FNORM)) MDEQC = 1
  655. MODE = MODE + MDEQC
  656. C
  657. C Check if solution to equality constraints satisfies inequality
  658. C constraints when there are no degrees of freedom left.
  659. C
  660. IF (KRANKE.EQ.N .AND. MG.GT.0) THEN
  661. XNORM = DASUM(N,X,1)
  662. MAPKE1 = MA + KRANKE + 1
  663. MEND = MA + KRANKE + MG
  664. DO 210 I = MAPKE1,MEND
  665. SIZE = DASUM(N,W(I,1),MDW)*XNORM + ABS(W(I,NP1))
  666. IF (W(I,NP1).GT.TAU*SIZE) THEN
  667. MODE = MODE + 2
  668. GO TO 290
  669. ENDIF
  670. 210 CONTINUE
  671. ENDIF
  672. ENDIF
  673. C
  674. C Replace diagonal terms of lower trapezoidal matrix.
  675. C
  676. IF (KRANKE.GT.0) THEN
  677. CALL DCOPY (KRANKE, WS(KRANKE+1), 1, W, MDW+1)
  678. C
  679. C Reapply transformation to put solution in original coordinates.
  680. C
  681. DO 220 I = KRANKE,1,-1
  682. CALL DH12 (2, I, I+1, N, W(I,1), MDW, WS(I), X, 1, 1, 1)
  683. 220 CONTINUE
  684. C
  685. C Compute covariance matrix of equality constrained problem.
  686. C
  687. IF (COV) THEN
  688. DO 270 J = MIN(KRANKE,N-1),1,-1
  689. RB = WS(J)*W(J,J)
  690. IF (RB.NE.0.D0) RB = 1.D0/RB
  691. JP1 = J + 1
  692. DO 230 I = JP1,N
  693. W(I,J) = RB*DDOT(N-J,W(I,JP1),MDW,W(J,JP1),MDW)
  694. 230 CONTINUE
  695. C
  696. GAM = 0.5D0*RB*DDOT(N-J,W(JP1,J),1,W(J,JP1),MDW)
  697. CALL DAXPY (N-J, GAM, W(J,JP1), MDW, W(JP1,J), 1)
  698. DO 250 I = JP1,N
  699. DO 240 K = I,N
  700. W(I,K) = W(I,K) + W(J,I)*W(K,J) + W(I,J)*W(J,K)
  701. W(K,I) = W(I,K)
  702. 240 CONTINUE
  703. 250 CONTINUE
  704. UJ = WS(J)
  705. VJ = GAM*UJ
  706. W(J,J) = UJ*VJ + UJ*VJ
  707. DO 260 I = JP1,N
  708. W(J,I) = UJ*W(I,J) + VJ*W(J,I)
  709. 260 CONTINUE
  710. CALL DCOPY (N-J, W(J, JP1), MDW, W(JP1,J), 1)
  711. 270 CONTINUE
  712. ENDIF
  713. ENDIF
  714. C
  715. C Apply the scaling to the covariance matrix.
  716. C
  717. IF (COV) THEN
  718. DO 280 I = 1,N
  719. CALL DSCAL (N, WS(I+N1-1), W(I,1), MDW)
  720. CALL DSCAL (N, WS(I+N1-1), W(1,I), 1)
  721. 280 CONTINUE
  722. ENDIF
  723. C
  724. C Rescale solution vector.
  725. C
  726. 290 IF (MODE.LE.1) THEN
  727. DO 300 J = 1,N
  728. X(J) = X(J)*WS(N1+J-1)
  729. 300 CONTINUE
  730. ENDIF
  731. C
  732. IP(1) = KRANKE
  733. IP(3) = IP(3) + 2*KRANKE + N
  734. RETURN
  735. END