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