1*DECK DWNLSM
2      SUBROUTINE DWNLSM (W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE,
3     +   IPIVOT, ITYPE, WD, H, SCALE, Z, TEMP, D)
4C***BEGIN PROLOGUE  DWNLSM
5C***SUBSIDIARY
6C***PURPOSE  Subsidiary to DWNNLS
7C***LIBRARY   SLATEC
8C***TYPE      DOUBLE PRECISION (WNLSM-S, DWNLSM-D)
9C***AUTHOR  Hanson, R. J., (SNLA)
10C           Haskell, K. H., (SNLA)
11C***DESCRIPTION
12C
13C     This is a companion subprogram to DWNNLS.
14C     The documentation for DWNNLS has complete usage instructions.
15C
16C     In addition to the parameters discussed in the prologue to
17C     subroutine DWNNLS, the following work arrays are used in
18C     subroutine DWNLSM  (they are passed through the calling
19C     sequence from DWNNLS for purposes of variable dimensioning).
20C     Their contents will in general be of no interest to the user.
21C
22C     Variables of type REAL are DOUBLE PRECISION.
23C
24C         IPIVOT(*)
25C            An array of length N.  Upon completion it contains the
26C         pivoting information for the cols of W(*,*).
27C
28C         ITYPE(*)
29C            An array of length M which is used to keep track
30C         of the classification of the equations.  ITYPE(I)=0
31C         denotes equation I as an equality constraint.
32C         ITYPE(I)=1 denotes equation I as a least squares
33C         equation.
34C
35C         WD(*)
36C            An array of length N.  Upon completion it contains the
37C         dual solution vector.
38C
39C         H(*)
40C            An array of length N.  Upon completion it contains the
41C         pivot scalars of the Householder transformations performed
42C         in the case KRANK.LT.L.
43C
44C         SCALE(*)
45C            An array of length M which is used by the subroutine
46C         to store the diagonal matrix of weights.
47C         These are used to apply the modified Givens
48C         transformations.
49C
50C         Z(*),TEMP(*)
51C            Working arrays of length N.
52C
53C         D(*)
54C            An array of length N that contains the
55C         column scaling for the matrix (E).
56C                                       (A)
57C
58C***SEE ALSO  DWNNLS
59C***ROUTINES CALLED  D1MACH, DASUM, DAXPY, DCOPY, DH12, DNRM2, DROTM,
60C                    DROTMG, DSCAL, DSWAP, DWNLIT, IDAMAX, XERMSG
61C***REVISION HISTORY  (YYMMDD)
62C   790701  DATE WRITTEN
63C   890531  Changed all specific intrinsics to generic.  (WRB)
64C   890618  Completely restructured and revised.  (WRB & RWC)
65C   891214  Prologue converted to Version 4.0 format.  (BAB)
66C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
67C   900328  Added TYPE section.  (WRB)
68C   900510  Fixed an error message.  (RWC)
69C   900604  DP version created from SP version.  (RWC)
70C   900911  Restriction on value of ALAMDA included.  (WRB)
71C***END PROLOGUE  DWNLSM
72      INTEGER IPIVOT(*), ITYPE(*), L, MA, MDW, MME, MODE, N
73      DOUBLE PRECISION D(*), H(*), PRGOPT(*), RNORM, SCALE(*), TEMP(*),
74     *   W(MDW,*), WD(*), X(*), Z(*)
75C
76      EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DH12, DNRM2, DROTM, DROTMG,
77     *   DSCAL, DSWAP, DWNLIT, IDAMAX, XERMSG
78      DOUBLE PRECISION D1MACH, DASUM, DNRM2
79      INTEGER IDAMAX
80C
81      DOUBLE PRECISION ALAMDA, ALPHA, ALSQ, AMAX, BLOWUP, BNORM,
82     *   DOPE(3), DRELPR, EANORM, FAC, SM, SPARAM(5), T, TAU, WMAX, Z2,
83     *   ZZ
84      INTEGER I, IDOPE(3), IMAX, ISOL, ITEMP, ITER, ITMAX, IWMAX, J,
85     *   JCON, JP, KEY, KRANK, L1, LAST, LINK, M, ME, NEXT, NIV, NLINK,
86     *   NOPT, NSOLN, NTIMES
87      LOGICAL DONE, FEASBL, FIRST, HITCON, POS
88C
89      SAVE DRELPR, FIRST
90      DATA FIRST /.TRUE./
91C***FIRST EXECUTABLE STATEMENT  DWNLSM
92C
93C     Initialize variables.
94C     DRELPR is the precision for the particular machine
95C     being used.  This logic avoids resetting it every entry.
96C
97      IF (FIRST) DRELPR = D1MACH(4)
98      FIRST = .FALSE.
99C
100C     Set the nominal tolerance used in the code.
101C
102      TAU = SQRT(DRELPR)
103C
104      M = MA + MME
105      ME = MME
106      MODE = 2
107C
108C     To process option vector
109C
110      FAC = 1.D-4
111C
112C     Set the nominal blow up factor used in the code.
113C
114      BLOWUP = TAU
115C
116C     The nominal column scaling used in the code is
117C     the identity scaling.
118C
119      CALL DCOPY (N, 1.D0, 0, D, 1)
120C
121C     Define bound for number of options to change.
122C
123      NOPT = 1000
124C
125C     Define bound for positive value of LINK.
126C
127      NLINK = 100000
128      NTIMES = 0
129      LAST = 1
130      LINK = PRGOPT(1)
131      IF (LINK.LE.0 .OR. LINK.GT.NLINK) THEN
132         CALL XERMSG ('SLATEC', 'DWNLSM',
133     +      'IN DWNNLS, THE OPTION VECTOR IS UNDEFINED', 3, 1)
134         RETURN
135      ENDIF
136C
137  100 IF (LINK.GT.1) THEN
138         NTIMES = NTIMES + 1
139         IF (NTIMES.GT.NOPT) THEN
140         CALL XERMSG ('SLATEC', 'DWNLSM',
141     +      'IN DWNNLS, THE LINKS IN THE OPTION VECTOR ARE CYCLING.',
142     +      3, 1)
143            RETURN
144         ENDIF
145C
146         KEY = PRGOPT(LAST+1)
147         IF (KEY.EQ.6 .AND. PRGOPT(LAST+2).NE.0.D0) THEN
148            DO 110 J = 1,N
149               T = DNRM2(M,W(1,J),1)
150               IF (T.NE.0.D0) T = 1.D0/T
151               D(J) = T
152  110       CONTINUE
153         ENDIF
154C
155         IF (KEY.EQ.7) CALL DCOPY (N, PRGOPT(LAST+2), 1, D, 1)
156         IF (KEY.EQ.8) TAU = MAX(DRELPR,PRGOPT(LAST+2))
157         IF (KEY.EQ.9) BLOWUP = MAX(DRELPR,PRGOPT(LAST+2))
158C
159         NEXT = PRGOPT(LINK)
160         IF (NEXT.LE.0 .OR. NEXT.GT.NLINK) THEN
161            CALL XERMSG ('SLATEC', 'DWNLSM',
162     +         'IN DWNNLS, THE OPTION VECTOR IS UNDEFINED', 3, 1)
163            RETURN
164         ENDIF
165C
166         LAST = LINK
167         LINK = NEXT
168         GO TO 100
169      ENDIF
170C
171      DO 120 J = 1,N
172         CALL DSCAL (M, D(J), W(1,J), 1)
173  120 CONTINUE
174C
175C     Process option vector
176C
177      DONE = .FALSE.
178      ITER = 0
179      ITMAX = 3*(N-L)
180      MODE = 0
181      NSOLN = L
182      L1 = MIN(M,L)
183C
184C     Compute scale factor to apply to equality constraint equations.
185C
186      DO 130 J = 1,N
187         WD(J) = DASUM(M,W(1,J),1)
188  130 CONTINUE
189C
190      IMAX = IDAMAX(N,WD,1)
191      EANORM = WD(IMAX)
192      BNORM = DASUM(M,W(1,N+1),1)
193      ALAMDA = EANORM/(DRELPR*FAC)
194C
195C     On machines, such as the VAXes using D floating, with a very
196C     limited exponent range for double precision values, the previously
197C     computed value of ALAMDA may cause an overflow condition.
198C     Therefore, this code further limits the value of ALAMDA.
199C
200      ALAMDA = MIN(ALAMDA,SQRT(D1MACH(2)))
201C
202C     Define scaling diagonal matrix for modified Givens usage and
203C     classify equation types.
204C
205      ALSQ = ALAMDA**2
206      DO 140 I = 1,M
207C
208C        When equation I is heavily weighted ITYPE(I)=0,
209C        else ITYPE(I)=1.
210C
211         IF (I.LE.ME) THEN
212            T = ALSQ
213            ITEMP = 0
214         ELSE
215            T = 1.D0
216            ITEMP = 1
217         ENDIF
218         SCALE(I) = T
219         ITYPE(I) = ITEMP
220  140 CONTINUE
221C
222C     Set the solution vector X(*) to zero and the column interchange
223C     matrix to the identity.
224C
225      CALL DCOPY (N, 0.D0, 0, X, 1)
226      DO 150 I = 1,N
227         IPIVOT(I) = I
228  150 CONTINUE
229C
230C     Perform initial triangularization in the submatrix
231C     corresponding to the unconstrained variables.
232C     Set first L components of dual vector to zero because
233C     these correspond to the unconstrained variables.
234C
235      CALL DCOPY (L, 0.D0, 0, WD, 1)
236C
237C     The arrays IDOPE(*) and DOPE(*) are used to pass
238C     information to DWNLIT().  This was done to avoid
239C     a long calling sequence or the use of COMMON.
240C
241      IDOPE(1) = ME
242      IDOPE(2) = NSOLN
243      IDOPE(3) = L1
244C
245      DOPE(1) = ALSQ
246      DOPE(2) = EANORM
247      DOPE(3) = TAU
248      CALL DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM,
249     +            IDOPE, DOPE, DONE)
250      ME    = IDOPE(1)
251      KRANK = IDOPE(2)
252      NIV   = IDOPE(3)
253C
254C     Perform WNNLS algorithm using the following steps.
255C
256C     Until(DONE)
257C        compute search direction and feasible point
258C        when (HITCON) add constraints
259C        else perform multiplier test and drop a constraint
260C        fin
261C     Compute-Final-Solution
262C
263C     To compute search direction and feasible point,
264C     solve the triangular system of currently non-active
265C     variables and store the solution in Z(*).
266C
267C     To solve system
268C     Copy right hand side into TEMP vector to use overwriting method.
269C
270  160 IF (DONE) GO TO 330
271      ISOL = L + 1
272      IF (NSOLN.GE.ISOL) THEN
273         CALL DCOPY (NIV, W(1,N+1), 1, TEMP, 1)
274         DO 170 J = NSOLN,ISOL,-1
275            IF (J.GT.KRANK) THEN
276               I = NIV - NSOLN + J
277            ELSE
278               I = J
279            ENDIF
280C
281            IF (J.GT.KRANK .AND. J.LE.L) THEN
282               Z(J) = 0.D0
283            ELSE
284               Z(J) = TEMP(I)/W(I,J)
285               CALL DAXPY (I-1, -Z(J), W(1,J), 1, TEMP, 1)
286            ENDIF
287  170    CONTINUE
288      ENDIF
289C
290C     Increment iteration counter and check against maximum number
291C     of iterations.
292C
293      ITER = ITER + 1
294      IF (ITER.GT.ITMAX) THEN
295         MODE = 1
296         DONE = .TRUE.
297      ENDIF
298C
299C     Check to see if any constraints have become active.
300C     If so, calculate an interpolation factor so that all
301C     active constraints are removed from the basis.
302C
303      ALPHA = 2.D0
304      HITCON = .FALSE.
305      DO 180 J = L+1,NSOLN
306         ZZ = Z(J)
307         IF (ZZ.LE.0.D0) THEN
308            T = X(J)/(X(J)-ZZ)
309            IF (T.LT.ALPHA) THEN
310               ALPHA = T
311               JCON = J
312            ENDIF
313            HITCON = .TRUE.
314         ENDIF
315  180 CONTINUE
316C
317C     Compute search direction and feasible point
318C
319      IF (HITCON) THEN
320C
321C        To add constraints, use computed ALPHA to interpolate between
322C        last feasible solution X(*) and current unconstrained (and
323C        infeasible) solution Z(*).
324C
325         DO 190 J = L+1,NSOLN
326            X(J) = X(J) + ALPHA*(Z(J)-X(J))
327  190    CONTINUE
328         FEASBL = .FALSE.
329C
330C        Remove column JCON and shift columns JCON+1 through N to the
331C        left.  Swap column JCON into the N th position.  This achieves
332C        upper Hessenberg form for the nonactive constraints and
333C        leaves an upper Hessenberg matrix to retriangularize.
334C
335  200    DO 210 I = 1,M
336            T = W(I,JCON)
337            CALL DCOPY (N-JCON, W(I, JCON+1), MDW, W(I, JCON), MDW)
338            W(I,N) = T
339  210    CONTINUE
340C
341C        Update permuted index vector to reflect this shift and swap.
342C
343         ITEMP = IPIVOT(JCON)
344         DO 220 I = JCON,N - 1
345            IPIVOT(I) = IPIVOT(I+1)
346  220    CONTINUE
347         IPIVOT(N) = ITEMP
348C
349C        Similarly permute X(*) vector.
350C
351         CALL DCOPY (N-JCON, X(JCON+1), 1, X(JCON), 1)
352         X(N) = 0.D0
353         NSOLN = NSOLN - 1
354         NIV = NIV - 1
355C
356C        Retriangularize upper Hessenberg matrix after adding
357C        constraints.
358C
359         I = KRANK + JCON - L
360         DO 230 J = JCON,NSOLN
361            IF (ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.0) THEN
362C
363C              Zero IP1 to I in column J
364C
365               IF (W(I+1,J).NE.0.D0) THEN
366                  CALL DROTMG (SCALE(I), SCALE(I+1), W(I,J), W(I+1,J),
367     +                         SPARAM)
368                  W(I+1,J) = 0.D0
369                  CALL DROTM (N+1-J, W(I,J+1), MDW, W(I+1,J+1), MDW,
370     +                        SPARAM)
371               ENDIF
372            ELSEIF (ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.1) THEN
373C
374C              Zero IP1 to I in column J
375C
376               IF (W(I+1,J).NE.0.D0) THEN
377                  CALL DROTMG (SCALE(I), SCALE(I+1), W(I,J), W(I+1,J),
378     +                         SPARAM)
379                  W(I+1,J) = 0.D0
380                  CALL DROTM (N+1-J, W(I,J+1), MDW, W(I+1,J+1), MDW,
381     +                        SPARAM)
382               ENDIF
383            ELSEIF (ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.0) THEN
384               CALL DSWAP (N+1, W(I,1), MDW, W(I+1,1), MDW)
385               CALL DSWAP (1, SCALE(I), 1, SCALE(I+1), 1)
386               ITEMP = ITYPE(I+1)
387               ITYPE(I+1) = ITYPE(I)
388               ITYPE(I) = ITEMP
389C
390C              Swapped row was formerly a pivot element, so it will
391C              be large enough to perform elimination.
392C              Zero IP1 to I in column J.
393C
394               IF (W(I+1,J).NE.0.D0) THEN
395                  CALL DROTMG (SCALE(I), SCALE(I+1), W(I,J), W(I+1,J),
396     +                         SPARAM)
397                  W(I+1,J) = 0.D0
398                  CALL DROTM (N+1-J, W(I,J+1), MDW, W(I+1,J+1), MDW,
399     +                        SPARAM)
400               ENDIF
401            ELSEIF (ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.1) THEN
402               IF (SCALE(I)*W(I,J)**2/ALSQ.GT.(TAU*EANORM)**2) THEN
403C
404C                 Zero IP1 to I in column J
405C
406                  IF (W(I+1,J).NE.0.D0) THEN
407                     CALL DROTMG (SCALE(I), SCALE(I+1), W(I,J),
408     +                            W(I+1,J), SPARAM)
409                     W(I+1,J) = 0.D0
410                     CALL DROTM (N+1-J, W(I,J+1), MDW, W(I+1,J+1), MDW,
411     +                           SPARAM)
412                  ENDIF
413               ELSE
414                  CALL DSWAP (N+1, W(I,1), MDW, W(I+1,1), MDW)
415                  CALL DSWAP (1, SCALE(I), 1, SCALE(I+1), 1)
416                  ITEMP = ITYPE(I+1)
417                  ITYPE(I+1) = ITYPE(I)
418                  ITYPE(I) = ITEMP
419                  W(I+1,J) = 0.D0
420               ENDIF
421            ENDIF
422            I = I + 1
423  230    CONTINUE
424C
425C        See if the remaining coefficients in the solution set are
426C        feasible.  They should be because of the way ALPHA was
427C        determined.  If any are infeasible, it is due to roundoff
428C        error.  Any that are non-positive will be set to zero and
429C        removed from the solution set.
430C
431         DO 240 JCON = L+1,NSOLN
432            IF (X(JCON).LE.0.D0) GO TO 250
433  240    CONTINUE
434         FEASBL = .TRUE.
435  250    IF (.NOT.FEASBL) GO TO 200
436      ELSE
437C
438C        To perform multiplier test and drop a constraint.
439C
440         CALL DCOPY (NSOLN, Z, 1, X, 1)
441         IF (NSOLN.LT.N) CALL DCOPY (N-NSOLN, 0.D0, 0, X(NSOLN+1), 1)
442C
443C        Reclassify least squares equations as equalities as necessary.
444C
445         I = NIV + 1
446  260    IF (I.LE.ME) THEN
447            IF (ITYPE(I).EQ.0) THEN
448               I = I + 1
449            ELSE
450               CALL DSWAP (N+1, W(I,1), MDW, W(ME,1), MDW)
451               CALL DSWAP (1, SCALE(I), 1, SCALE(ME), 1)
452               ITEMP = ITYPE(I)
453               ITYPE(I) = ITYPE(ME)
454               ITYPE(ME) = ITEMP
455               ME = ME - 1
456            ENDIF
457            GO TO 260
458         ENDIF
459C
460C        Form inner product vector WD(*) of dual coefficients.
461C
462         DO 280 J = NSOLN+1,N
463            SM = 0.D0
464            DO 270 I = NSOLN+1,M
465               SM = SM + SCALE(I)*W(I,J)*W(I,N+1)
466  270       CONTINUE
467            WD(J) = SM
468  280    CONTINUE
469C
470C        Find J such that WD(J)=WMAX is maximum.  This determines
471C        that the incoming column J will reduce the residual vector
472C        and be positive.
473C
474  290    WMAX = 0.D0
475         IWMAX = NSOLN + 1
476         DO 300 J = NSOLN+1,N
477            IF (WD(J).GT.WMAX) THEN
478               WMAX = WD(J)
479               IWMAX = J
480            ENDIF
481  300    CONTINUE
482         IF (WMAX.LE.0.D0) GO TO 330
483C
484C        Set dual coefficients to zero for incoming column.
485C
486         WD(IWMAX) = 0.D0
487C
488C        WMAX .GT. 0.D0, so okay to move column IWMAX to solution set.
489C        Perform transformation to retriangularize, and test for near
490C        linear dependence.
491C
492C        Swap column IWMAX into NSOLN-th position to maintain upper
493C        Hessenberg form of adjacent columns, and add new column to
494C        triangular decomposition.
495C
496         NSOLN = NSOLN + 1
497         NIV = NIV + 1
498         IF (NSOLN.NE.IWMAX) THEN
499            CALL DSWAP (M, W(1,NSOLN), 1, W(1,IWMAX), 1)
500            WD(IWMAX) = WD(NSOLN)
501            WD(NSOLN) = 0.D0
502            ITEMP = IPIVOT(NSOLN)
503            IPIVOT(NSOLN) = IPIVOT(IWMAX)
504            IPIVOT(IWMAX) = ITEMP
505         ENDIF
506C
507C        Reduce column NSOLN so that the matrix of nonactive constraints
508C        variables is triangular.
509C
510         DO 320 J = M,NIV+1,-1
511            JP = J - 1
512C
513C           When operating near the ME line, test to see if the pivot
514C           element is near zero.  If so, use the largest element above
515C           it as the pivot.  This is to maintain the sharp interface
516C           between weighted and non-weighted rows in all cases.
517C
518            IF (J.EQ.ME+1) THEN
519               IMAX = ME
520               AMAX = SCALE(ME)*W(ME,NSOLN)**2
521               DO 310 JP = J - 1,NIV,-1
522                  T = SCALE(JP)*W(JP,NSOLN)**2
523                  IF (T.GT.AMAX) THEN
524                     IMAX = JP
525                     AMAX = T
526                  ENDIF
527  310          CONTINUE
528               JP = IMAX
529            ENDIF
530C
531            IF (W(J,NSOLN).NE.0.D0) THEN
532               CALL DROTMG (SCALE(JP), SCALE(J), W(JP,NSOLN),
533     +                      W(J,NSOLN), SPARAM)
534               W(J,NSOLN) = 0.D0
535               CALL DROTM (N+1-NSOLN, W(JP,NSOLN+1), MDW, W(J,NSOLN+1),
536     +                     MDW, SPARAM)
537            ENDIF
538  320    CONTINUE
539C
540C        Solve for Z(NSOLN)=proposed new value for X(NSOLN).  Test if
541C        this is nonpositive or too large.  If this was true or if the
542C        pivot term was zero, reject the column as dependent.
543C
544         IF (W(NIV,NSOLN).NE.0.D0) THEN
545            ISOL = NIV
546            Z2 = W(ISOL,N+1)/W(ISOL,NSOLN)
547            Z(NSOLN) = Z2
548            POS = Z2 .GT. 0.D0
549            IF (Z2*EANORM.GE.BNORM .AND. POS) THEN
550               POS = .NOT. (BLOWUP*Z2*EANORM.GE.BNORM)
551            ENDIF
552C
553C           Try to add row ME+1 as an additional equality constraint.
554C           Check size of proposed new solution component.
555C           Reject it if it is too large.
556C
557         ELSEIF (NIV.LE.ME .AND. W(ME+1,NSOLN).NE.0.D0) THEN
558            ISOL = ME + 1
559            IF (POS) THEN
560C
561C              Swap rows ME+1 and NIV, and scale factors for these rows.
562C
563               CALL DSWAP (N+1, W(ME+1,1), MDW, W(NIV,1), MDW)
564               CALL DSWAP (1, SCALE(ME+1), 1, SCALE(NIV), 1)
565               ITEMP = ITYPE(ME+1)
566               ITYPE(ME+1) = ITYPE(NIV)
567               ITYPE(NIV) = ITEMP
568               ME = ME + 1
569            ENDIF
570         ELSE
571            POS = .FALSE.
572         ENDIF
573C
574         IF (.NOT.POS) THEN
575            NSOLN = NSOLN - 1
576            NIV = NIV - 1
577         ENDIF
578         IF (.NOT.(POS.OR.DONE)) GO TO 290
579      ENDIF
580      GO TO 160
581C
582C     Else perform multiplier test and drop a constraint.  To compute
583C     final solution.  Solve system, store results in X(*).
584C
585C     Copy right hand side into TEMP vector to use overwriting method.
586C
587  330 ISOL = 1
588      IF (NSOLN.GE.ISOL) THEN
589         CALL DCOPY (NIV, W(1,N+1), 1, TEMP, 1)
590         DO 340 J = NSOLN,ISOL,-1
591            IF (J.GT.KRANK) THEN
592               I = NIV - NSOLN + J
593            ELSE
594               I = J
595            ENDIF
596C
597            IF (J.GT.KRANK .AND. J.LE.L) THEN
598               Z(J) = 0.D0
599            ELSE
600               Z(J) = TEMP(I)/W(I,J)
601               CALL DAXPY (I-1, -Z(J), W(1,J), 1, TEMP, 1)
602            ENDIF
603  340    CONTINUE
604      ENDIF
605C
606C     Solve system.
607C
608      CALL DCOPY (NSOLN, Z, 1, X, 1)
609C
610C     Apply Householder transformations to X(*) if KRANK.LT.L
611C
612      IF (KRANK.LT.L) THEN
613         DO 350 I = 1,KRANK
614            CALL DH12 (2, I, KRANK+1, L, W(I,1), MDW, H(I), X, 1, 1, 1)
615  350    CONTINUE
616      ENDIF
617C
618C     Fill in trailing zeroes for constrained variables not in solution.
619C
620      IF (NSOLN.LT.N) CALL DCOPY (N-NSOLN, 0.D0, 0, X(NSOLN+1), 1)
621C
622C     Permute solution vector to natural order.
623C
624      DO 380 I = 1,N
625         J = I
626  360    IF (IPIVOT(J).EQ.I) GO TO 370
627         J = J + 1
628         GO TO 360
629C
630  370    IPIVOT(J) = IPIVOT(I)
631         IPIVOT(I) = J
632         CALL DSWAP (1, X(J), 1, X(I), 1)
633  380 CONTINUE
634C
635C     Rescale the solution using the column scaling.
636C
637      DO 390 J = 1,N
638         X(J) = X(J)*D(J)
639  390 CONTINUE
640C
641      DO 400 I = NSOLN+1,M
642         T = W(I,N+1)
643         IF (I.LE.ME) T = T/ALAMDA
644         T = (SCALE(I)*T)*T
645         RNORM = RNORM + T
646  400 CONTINUE
647C
648      RNORM = SQRT(RNORM)
649      RETURN
650      END
651