SUBROUTINE LSQSL2 (NDIM,A,D,W,B,X,IRANK,IN,ITMAX,IT,IEQ,ENORM,EPS1 2LSQLS2.38
1,NHDIM,H,AA,R,S) LSQLS2.39
C LSQLS2.40
C RH141293.94
IMPLICIT NONE RH141293.95
C RH141293.96
C LSQLS2.41
LOGICAL ERM LSQLS2.42
INTEGER D,W LSQLS2.43
& ,I ! Loop index RH141293.97
& ,IEQ ! Column scaling indicator RH141293.98
& ,II ! IP - 1 RH141293.99
& ,IN ! Least squares option no. RH141293.100
& ,IP ! Pivot row RH141293.101
& ,IPM1 ! Pivot row - 1 RH141293.102
& ,IPP1 ! Pivot row + 1 RH141293.103
& ,IRANK ! RH141293.104
& ,IRM1 ! IRANK - 1 RH141293.105
& ,IRP1 ! IRANK + 1 RH141293.106
& ,ISW ! RH141293.107
& ,IT ! Iteration counter RH141293.108
& ,ITMAX ! Max no of iterations RH141293.109
& ,J ! Loop index RH141293.110
& ,J1,J2,J3,J4,J5,J6,J7,J8,J9 ! RH141293.111
& ,K1 ! RH141293.112
& ,K ! Row index for A RH141293.113
& ,L ! Loop index RH141293.114
& ,LM ! RH141293.115
& ,N ! Column index for A RH141293.116
& ,NS ! RH141293.117
& ,N1,N2,N3,N4,N5,N6,N7,N8 ! RH141293.118
& ,NDIM ! Dimension for arrays RH141293.119
& ,NHDIM ! Dimension for arrays RH141293.120
& ,M ! RH141293.121
C LSQLS2.44
C IN=1 FOR FIRST ENTRY. LSQLS2.45
C A IS DECOMPOSED AND SAVED. AX-B IS SOLVED. LSQLS2.46
C IN = 2 FOR SUBSEQUENT ENTRIES WITH A NEW VECTOR B. LSQLS2.47
C IN=3 TO RESTORE A FROM THE PREVIOUS ENTRY. LSQLS2.48
C IN=4 TO CONTINUE THE ITERATIVE IMPROVEMENT FOR THIS SYSTEM. LSQLS2.49
C IN = 5 TO CALCULATE SOLUTIONS TO AX=0, THEN STORE IN THE ARRAY H. LSQLS2.50
C IN = 6 DO NOT STORE A IN AA. OBTAIN T = QAS, WHERE T IS LSQLS2.51
C MIN(K,N) X MIN(K,N) AND UPPER TRIANGULAR. NOW RETURN.DO NOT OBTAIN LSQLS2.52
C A SOLUTION. LSQLS2.53
C NO SCALING OR COLUMN INTERCHANGES ARE PERFORMED. LSQLS2.54
C IN = 7 SAME AS WITH IN = 6 EXCEPT THAT SOLN. OF MIN. LENGTH LSQLS2.55
C IS PLACED INTO X. NO ITERATIVE REFINEMENT. NOW RETURN. LSQLS2.56
C COLUMN INTERCHANGES ARE PERFORMED. NO SCALING IS PERFORMED. LSQLS2.57
C IN = 8 SET ADDRESSES. NOW RETURN. LSQLS2.58
C LSQLS2.59
C OPTIONS FOR COMPUTING A MATRIX PRODUCT Y*H OR H*Y ARE LSQLS2.60
C AVAILABLE WITH THE USE OF THE ENTRY POINTS MYH AND MHY. LSQLS2.61
C USE OF THESE OPTIONS IN THESE ENTRY POINTS ALLOW A GREAT SAVING IN LSQLS2.62
C STORAGE REQUIRED. LSQLS2.63
C LSQLS2.64
C LSQLS2.65
REAL A(NDIM,NDIM),AA(D,W),H(NHDIM,NHDIM) @DYALLOC.4435
REAL B(NDIM),S(4*NDIM), X(NHDIM),R(4*NDIM) @DYALLOC.4436
& ,A1 ! RH141293.122
& ,A2 ! RH141293.123
& ,AJ ! Jth row of A RH141293.124
& ,AM ! RH141293.125
& ,AZ ! RH141293.126
& ,BP ! 2/length of U**2 RH141293.127
& ,DP ! RH141293.128
& ,EN1 ! RH141293.129
& ,ENM1 ! RH141293.130
& ,ENORM ! RH141293.131
& ,EPS1 ! Rank control variable RH141293.132
& ,EPS2 ! Rank comtrol variable RH141293.133
& ,UP ! RH141293.134
& ,SJ ! RH141293.135
& ,SP ! RH141293.136
& ,TOP ! RH141293.137
& ,TOP1 ! RH141293.138
& ,TOP2 ! RH141293.139
SAVE IRP1 OSI1F402.23
C D = DEPTH OF MATRIX. LSQLS2.67
C W = WIDTH OF MATRIX. LSQLS2.68
K=D LSQLS2.69
N=W LSQLS2.70
ERM=.TRUE. LSQLS2.71
C LSQLS2.72
C IF IT=0 ON ENTRY, THE POSSIBLE ERROR MESSAGE WILL BE SUPPRESSED. LSQLS2.73
C LSQLS2.74
IF (IT.EQ.0) ERM=.FALSE. LSQLS2.75
C LSQLS2.76
C IEQ = 2 IF COLUMN SCALING BY LEAST MAX. COLUMN LENGTH IS LSQLS2.77
C TO BE PERFORMED. LSQLS2.78
C LSQLS2.79
C IEQ = 1 IF SCALING OF ALL COMPONENTS IS TO BE DONE WITH LSQLS2.80
C THE SCALAR MAX(ABS(AIJ))/K*N. LSQLS2.81
C LSQLS2.82
C IEQ = 3 IF COLUMN SCALING AS WITH IN =2 WILL BE RETAINED IN LSQLS2.83
C RANK DEFICIENT CASES. LSQLS2.84
C LSQLS2.85
C THE ARRAY S MUST CONTAIN AT LEAST MAX(K,N) + 4N + 4MIN(K,N) CELLS LSQLS2.86
C THE ARRAY R MUST CONTAIN K+4N S.P. CELLS. LSQLS2.87
C LSQLS2.88
DATA EPS2/1.E-16/ LSQLS2.89
C THE LAST CARD CONTROLS DESIRED RELATIVE ACCURACY. LSQLS2.90
C EPS1 CONTROLS (EPS) RANK. LSQLS2.91
C LSQLS2.92
ISW=1 LSQLS2.93
L=MIN0(K,N) LSQLS2.94
M=MAX0(K,N) LSQLS2.95
J1=M LSQLS2.96
J2=N+J1 LSQLS2.97
J3=J2+N LSQLS2.98
J4=J3+L LSQLS2.99
J5=J4+L LSQLS2.100
J6=J5+L LSQLS2.101
J7=J6+L LSQLS2.102
J8=J7+N LSQLS2.103
J9=J8+N LSQLS2.104
LM=L LSQLS2.105
IF (IRANK.GE.1.AND.IRANK.LE.L) LM=IRANK LSQLS2.106
IF (IN.EQ.6) LM=L LSQLS2.107
IF (IN.EQ.8) RETURN LSQLS2.108
C LSQLS2.109
C RETURN AFTER SETTING ADDRESSES WHEN IN=8. LSQLS2.110
C LSQLS2.111
GO TO (10,360,810,390,830,10,10), IN LSQLS2.112
C LSQLS2.113
C EQUILIBRATE COLUMNS OF A (1)-(2). LSQLS2.114
C LSQLS2.115
C (1) LSQLS2.116
C LSQLS2.117
10 CONTINUE LSQLS2.118
C LSQLS2.119
C SAVE DATA WHEN IN = 1. LSQLS2.120
C LSQLS2.121
IF (IN.GT.5) GO TO 30 LSQLS2.122
DO 20 J=1,N LSQLS2.123
DO 20 I=1,K LSQLS2.124
20 AA(I,J)=A(I,J) LSQLS2.125
30 CONTINUE LSQLS2.126
IF (IEQ.EQ.1) GO TO 60 LSQLS2.127
DO 50 J=1,N LSQLS2.128
AM=0.E0 LSQLS2.129
DO 40 I=1,K LSQLS2.130
40 AM=MAX(AM,ABS(A(I,J))) LSQLS2.131
C LSQLS2.132
C S(M+N+1)-S(M+2N) CONTAINS SCALING FOR OUTPUT VARIABLES. LSQLS2.133
C LSQLS2.134
N2=J2+J LSQLS2.135
IF (IN.EQ.6) AM=1.E0 LSQLS2.136
S(N2)=1.E0/AM LSQLS2.137
DO 50 I=1,K LSQLS2.138
50 A(I,J)=A(I,J)*S(N2) LSQLS2.139
GO TO 100 LSQLS2.140
60 AM=0.E0 LSQLS2.141
DO 70 J=1,N LSQLS2.142
DO 70 I=1,K LSQLS2.143
70 AM=MAX(AM,ABS(A(I,J))) LSQLS2.144
AM=AM/FLOAT(K*N) LSQLS2.145
IF (IN.EQ.6) AM=1.E0 LSQLS2.146
DO 80 J=1,N LSQLS2.147
N2=J2+J LSQLS2.148
80 S(N2)=1.E0/AM LSQLS2.149
DO 90 J=1,N LSQLS2.150
N2=J2+J LSQLS2.151
DO 90 I=1,K LSQLS2.152
90 A(I,J)=A(I,J)*S(N2) LSQLS2.153
C COMPUTE COLUMN LENGTHS WITH D.P. SUMS FINALLY ROUNDED TO S.P. LSQLS2.154
C LSQLS2.155
C (2) LSQLS2.156
C LSQLS2.157
100 DO 110 J=1,N LSQLS2.158
N7=J7+J LSQLS2.159
N2=J2+J LSQLS2.160
110 S(N7)=S(N2) LSQLS2.161
C LSQLS2.162
C S(M+1)-S(M+ N) CONTAINS VARIABLE PERMUTATIONS. LSQLS2.163
C LSQLS2.164
C SET PERMUTATION TO IDENTITY. LSQLS2.165
C LSQLS2.166
DO 120 J=1,N LSQLS2.167
N1=J1+J LSQLS2.168
120 S(N1)=J LSQLS2.169
C LSQLS2.170
C BEGIN ELIMINATION ON THE MATRIX A WITH ORTHOGONAL MATRICES . LSQLS2.171
C LSQLS2.172
C IP=PIVOT ROW LSQLS2.173
C LSQLS2.174
DO 250 IP=1,LM LSQLS2.175
C LSQLS2.176
C LSQLS2.177
*IF -DEF,T3E OSI1F402.25
DP=0.D0 LSQLS2.178
*ELSE OSI1F402.26
DP=0.0 OSI1F402.27
*ENDIF OSI1F402.28
LM=IP LSQLS2.179
DO 140 J=IP,N LSQLS2.180
*IF -DEF,T3E OSI1F402.29
SJ=0.D0 LSQLS2.181
*ELSE OSI1F402.30
SJ=0.0 OSI1F402.31
*ENDIF OSI1F402.32
DO 130 I=IP,K LSQLS2.182
SJ=SJ+A(I,J)**2 LSQLS2.183
130 CONTINUE LSQLS2.184
IF (DP.GT.SJ) GO TO 140 LSQLS2.185
DP=SJ LSQLS2.186
LM=J LSQLS2.187
IF (IN.EQ.6) GO TO 160 LSQLS2.188
140 CONTINUE LSQLS2.189
C LSQLS2.190
C MAXIMIZE (SIGMA)**2 BY COLUMN INTERCHANGE. LSQLS2.191
C LSQLS2.192
C SUPRESS COLUMN INTERCHANGES WHEN IN=6. LSQLS2.193
C LSQLS2.194
C LSQLS2.195
C EXCHANGE COLUMNS IF NECESSARY. LSQLS2.196
C LSQLS2.197
IF (LM.EQ.IP) GO TO 160 LSQLS2.198
DO 150 I=1,K LSQLS2.199
A1=A(I,IP) LSQLS2.200
A(I,IP)=A(I,LM) LSQLS2.201
150 A(I,LM)=A1 LSQLS2.202
C LSQLS2.203
C RECORD PERMUTATION AND EXCHANGE SQUARES OF COLUMN LENGTHS. LSQLS2.204
C LSQLS2.205
N1=J1+LM LSQLS2.206
A1=S(N1) LSQLS2.207
N2=J1+IP LSQLS2.208
S(N1)=S(N2) LSQLS2.209
S(N2)=A1 LSQLS2.210
N7=J7+LM LSQLS2.211
N8=J7+IP LSQLS2.212
A1=S(N7) LSQLS2.213
S(N7)=S(N8) LSQLS2.214
S(N8)=A1 LSQLS2.215
160 IF (IP.EQ.1) GO TO 180 LSQLS2.216
A1=0.E0 LSQLS2.217
IPM1=IP-1 LSQLS2.218
DO 170 I=1,IPM1 LSQLS2.219
A1=A1+A(I,IP)**2 LSQLS2.220
170 CONTINUE LSQLS2.221
IF (A1.GT.0.E0) GO TO 190 LSQLS2.222
*IF -DEF,T3E OSI1F402.33
180 IF (DP.GT.0.D0) GO TO 200 LSQLS2.223
*ELSE OSI1F402.34
180 IF (DP.GT.0.0) GO TO 200 OSI1F402.35
*ENDIF OSI1F402.36
C LSQLS2.224
C TEST FOR RANK DEFICIENCY. LSQLS2.225
C LSQLS2.226
190 IF (SQRT(DP/A1).GT.EPS1) GO TO 200 LSQLS2.227
IF (IN.EQ.6) GO TO 200 LSQLS2.228
II=IP-1 LSQLS2.229
IF (ERM) WRITE (6,1140) IRANK,EPS1,II,II LSQLS2.230
IRANK=IP-1 LSQLS2.231
ERM=.FALSE. LSQLS2.232
GO TO 260 LSQLS2.233
C LSQLS2.234
C (EPS1) RANK IS DEFICIENT. LSQLS2.235
C LSQLS2.236
200 SP=SQRT(DP) LSQLS2.237
C LSQLS2.238
C BEGIN FRONT ELIMINATION ON COLUMN IP. LSQLS2.239
C LSQLS2.240
C SP=SQROOT(SIGMA**2). LSQLS2.241
C LSQLS2.242
*IF -DEF,T3E OSI1F402.37
BP=1.D0/(DP+SP*ABS(A(IP,IP))) LSQLS2.243
*ELSE OSI1F402.38
BP=1.0/(DP+SP*ABS(A(IP,IP))) OSI1F402.39
*ENDIF OSI1F402.40
C LSQLS2.244
C STORE BETA IN S(3N+1)-S(3N+L). LSQLS2.245
C LSQLS2.246
*IF -DEF,T3E OSI1F402.41
IF (IP.EQ.K) BP=0.D0 LSQLS2.247
*ELSE OSI1F402.42
IF (IP.EQ.K) BP=0.0 OSI1F402.43
*ENDIF OSI1F402.44
N3=K+2*N+IP LSQLS2.248
R(N3)=BP LSQLS2.249
UP=SIGN(SP+ABS(A(IP,IP)),A(IP,IP)) LSQLS2.250
IF (IP.GE.K) GO TO 250 LSQLS2.251
IPP1=IP+1 LSQLS2.252
IF (IP.GE.N) GO TO 240 LSQLS2.253
DO 230 J=IPP1,N LSQLS2.254
*IF -DEF,T3E OSI1F402.45
SJ=0.D0 LSQLS2.255
*ELSE OSI1F402.46
SJ=0.0 OSI1F402.47
*ENDIF OSI1F402.48
DO 210 I=IPP1,K LSQLS2.256
210 SJ=SJ+A(I,J)*A(I,IP) LSQLS2.257
SJ=SJ+UP*A(IP,J) LSQLS2.258
SJ=BP*SJ LSQLS2.259
C LSQLS2.260
C SJ=YJ NOW LSQLS2.261
C LSQLS2.262
DO 220 I=IPP1,K LSQLS2.263
220 A(I,J)=A(I,J)-A(I,IP)*SJ LSQLS2.264
230 A(IP,J)=A(IP,J)-SJ*UP LSQLS2.265
240 A(IP,IP)=-SIGN(SP,A(IP,IP)) LSQLS2.266
C LSQLS2.267
N4=K+3*N+IP LSQLS2.268
R(N4)=UP LSQLS2.269
250 CONTINUE LSQLS2.270
IRANK=LM LSQLS2.271
260 IRP1=IRANK+1 LSQLS2.272
IRM1=IRANK-1 LSQLS2.273
IF (IRANK.EQ.0.OR.IRANK.EQ.N) GO TO 360 LSQLS2.274
IF (IEQ.EQ.3) GO TO 290 LSQLS2.275
C LSQLS2.276
C BEGIN BACK PROCESSING FOR RANK DEFICIENCY CASE LSQLS2.277
C IF IRANK IS LESS THAN N. LSQLS2.278
C LSQLS2.279
DO 280 J=1,N LSQLS2.280
N2=J2+J LSQLS2.281
N7=J7+J LSQLS2.282
L=MIN0(J,IRANK) LSQLS2.283
C LSQLS2.284
C UNSCALE COLUMNS FOR RANK DEFICIENT MATRICES WHEN IEQ.NE.3. LSQLS2.285
C LSQLS2.286
DO 270 I=1,L LSQLS2.287
270 A(I,J)=A(I,J)/S(N7) LSQLS2.288
S(N7)=1.E0 LSQLS2.289
280 S(N2)=1.E0 LSQLS2.290
290 IP=IRANK LSQLS2.291
*IF -DEF,T3E OSI1F402.49
300 SJ=0.D0 LSQLS2.292
*ELSE OSI1F402.50
300 SJ=0.0 OSI1F402.51
*ENDIF OSI1F402.52
DO 310 J=IRP1,N LSQLS2.293
SJ=SJ+A(IP,J)**2 LSQLS2.294
310 CONTINUE LSQLS2.295
SJ=SJ+A(IP,IP)**2 LSQLS2.296
AJ=SQRT(SJ) LSQLS2.297
UP=SIGN(AJ+ABS(A(IP,IP)),A(IP,IP)) LSQLS2.298
C LSQLS2.299
C IP TH ELEMENT OF U VECTOR CALCULATED. LSQLS2.300
C LSQLS2.301
*IF -DEF,T3E OSI1F402.53
BP=1.D0/(SJ+ABS(A(IP,IP))*AJ) LSQLS2.302
*ELSE OSI1F402.54
BP=1.0/(SJ+ABS(A(IP,IP))*AJ) OSI1F402.55
*ENDIF OSI1F402.56
C LSQLS2.303
C BP = 2/LENGTH OF U SQUARED. LSQLS2.304
C LSQLS2.305
IPM1=IP-1 LSQLS2.306
IF (IPM1.LE.0) GO TO 340 LSQLS2.307
DO 330 I=1,IPM1 LSQLS2.308
DP=A(I,IP)*UP LSQLS2.309
DO 320 J=IRP1,N LSQLS2.310
DP=DP+A(I,J)*A(IP,J) LSQLS2.311
320 CONTINUE LSQLS2.312
DP=DP/(SJ+ABS(A(IP,IP))*AJ) LSQLS2.313
C LSQLS2.314
C CALC. (AJ,U), WHERE AJ=JTH ROW OF A LSQLS2.315
C LSQLS2.316
A(I,IP)=A(I,IP)-UP*DP LSQLS2.317
C LSQLS2.318
C MODIFY ARRAY A. LSQLS2.319
C LSQLS2.320
DO 330 J=IRP1,N LSQLS2.321
330 A(I,J)=A(I,J)-A(IP,J)*DP LSQLS2.322
340 A(IP,IP)=-SIGN(AJ,A(IP,IP)) LSQLS2.323
C LSQLS2.324
C CALC. MODIFIED PIVOT. LSQLS2.325
C LSQLS2.326
C LSQLS2.327
C SAVE BETA AND IP TH ELEMENT OF U VECTOR IN R ARRAY. LSQLS2.328
C LSQLS2.329
N6=K+IP LSQLS2.330
N7=K+N+IP LSQLS2.331
R(N6)=BP LSQLS2.332
R(N7)=UP LSQLS2.333
C LSQLS2.334
C TEST FOR END OF BACK PROCESSING. LSQLS2.335
C LSQLS2.336
IF (IP-1) 360,360,350 LSQLS2.337
350 IP=IP-1 LSQLS2.338
GO TO 300 LSQLS2.339
360 IF (IN.EQ.6) RETURN LSQLS2.340
DO 370 J=1,K LSQLS2.341
370 R(J)=B(J) LSQLS2.342
IT=0 LSQLS2.343
C LSQLS2.344
C SET INITIAL X VECTOR TO ZERO. LSQLS2.345
C LSQLS2.346
DO 380 J=1,N LSQLS2.347
*IF -DEF,T3E OSI1F402.57
380 X(J)=0.D0 LSQLS2.348
*ELSE OSI1F402.58
380 X(J)=0.0 OSI1F402.59
*ENDIF OSI1F402.60
IF (IRANK.EQ.0) GO TO 690 LSQLS2.349
C LSQLS2.350
C APPLY Q TO RT. HAND SIDE. LSQLS2.351
C LSQLS2.352
390 DO 430 IP=1,IRANK LSQLS2.353
N4=K+3*N+IP LSQLS2.354
SJ=R(N4)*R(IP) LSQLS2.355
IPP1=IP+1 LSQLS2.356
IF (IPP1.GT.K) GO TO 410 LSQLS2.357
DO 400 I=IPP1,K LSQLS2.358
400 SJ=SJ+A(I,IP)*R(I) LSQLS2.359
410 N3=K+2*N+IP LSQLS2.360
BP=R(N3) LSQLS2.361
IF (IPP1.GT.K) GO TO 430 LSQLS2.362
DO 420 I=IPP1,K LSQLS2.363
420 R(I)=R(I)-BP*A(I,IP)*SJ LSQLS2.364
430 R(IP)=R(IP)-BP*R(N4)*SJ LSQLS2.365
DO 440 J=1,IRANK LSQLS2.366
440 S(J)=R(J) LSQLS2.367
ENORM=0.E0 LSQLS2.368
IF (IRP1.GT.K) GO TO 510 LSQLS2.369
DO 450 J=IRP1,K LSQLS2.370
450 ENORM=ENORM+R(J)**2 LSQLS2.371
ENORM=SQRT(ENORM) LSQLS2.372
GO TO 510 LSQLS2.373
460 DO 480 J=1,N LSQLS2.374
*IF -DEF,T3E OSI1F402.61
SJ=0.D0 LSQLS2.375
*ELSE OSI1F402.62
SJ=0.0 OSI1F402.63
*ENDIF OSI1F402.64
N1=J1+J LSQLS2.376
IP=S(N1) LSQLS2.377
DO 470 I=1,K LSQLS2.378
470 SJ=SJ+R(I)*AA(I,IP) LSQLS2.379
C LSQLS2.380
C APPLY AT TO RT. HAND SIDE. LSQLS2.381
C APPLY SCALING. LSQLS2.382
C LSQLS2.383
N7=J2+IP LSQLS2.384
N1=K+N+J LSQLS2.385
480 R(N1)=SJ*S(N7) LSQLS2.386
N1=K+N LSQLS2.387
S(1)=R(N1+1)/A(1,1) LSQLS2.388
IF (N.EQ.1) GO TO 510 LSQLS2.389
DO 500 J=2,N LSQLS2.390
N1=J-1 LSQLS2.391
*IF -DEF,T3E OSI1F402.65
SJ=0.D0 LSQLS2.392
*ELSE OSI1F402.66
SJ=0.0 OSI1F402.67
*ENDIF OSI1F402.68
DO 490 I=1,N1 LSQLS2.393
490 SJ=SJ+A(I,J)*S(I) LSQLS2.394
N2=K+J+N LSQLS2.395
500 S(J)=(R(N2)-SJ)/A(J,J) LSQLS2.396
C LSQLS2.397
C ENTRY TO CONTINUE ITERATING. SOLVES TZ = C = 1ST IRANK LSQLS2.398
C COMPONENTS OF QB . LSQLS2.399
C LSQLS2.400
510 S(IRANK)=S(IRANK)/A(IRANK,IRANK) LSQLS2.401
IF (IRM1.EQ.0) GO TO 540 LSQLS2.402
DO 530 J=1,IRM1 LSQLS2.403
N1=IRANK-J LSQLS2.404
N2=N1+1 LSQLS2.405
SJ=0. LSQLS2.406
DO 520 I=N2,IRANK LSQLS2.407
520 SJ=SJ+A(N1,I)*S(I) LSQLS2.408
530 S(N1)=(S(N1)-SJ)/A(N1,N1) LSQLS2.409
C LSQLS2.410
C Z CALCULATED. COMPUTE X = SZ. LSQLS2.411
C LSQLS2.412
540 IF (IRANK.EQ.N) GO TO 590 LSQLS2.413
DO 550 J=IRP1,N LSQLS2.414
550 S(J)=0.E0 LSQLS2.415
DO 580 I=1,IRANK LSQLS2.416
N7=K+N+I LSQLS2.417
SJ=R(N7)*S(I) LSQLS2.418
DO 560 J=IRP1,N LSQLS2.419
SJ=SJ+A(I,J)*S(J) LSQLS2.420
560 CONTINUE LSQLS2.421
N6=K+I LSQLS2.422
DO 570 J=IRP1,N LSQLS2.423
570 S(J)=S(J)-A(I,J)*R(N6)*SJ LSQLS2.424
580 S(I)=S(I)-R(N6)*R(N7)*SJ LSQLS2.425
C LSQLS2.426
C INCREMENT FOR X OF MINIMAL LENGTH CALCULATED. LSQLS2.427
C LSQLS2.428
590 DO 600 I=1,N LSQLS2.429
600 X(I)=X(I)+S(I) LSQLS2.430
IF (IN.EQ.7) GO TO 750 LSQLS2.431
C LSQLS2.432
C CALC. SUP NORM OF INCREMENT AND RESIDUALS LSQLS2.433
C LSQLS2.434
TOP1=0.E0 LSQLS2.435
DO 610 J=1,N LSQLS2.436
N2=J7+J LSQLS2.437
610 TOP1=MAX(TOP1,ABS(S(J))*S(N2)) LSQLS2.438
DO 630 I=1,K LSQLS2.439
*IF -DEF,T3E OSI1F402.69
SJ=0.D0 LSQLS2.440
*ELSE OSI1F402.70
SJ=0.0 OSI1F402.71
*ENDIF OSI1F402.72
DO 620 J=1,N LSQLS2.441
N1=J1+J LSQLS2.442
IP=S(N1) LSQLS2.443
N7=J2+IP LSQLS2.444
620 SJ=SJ+AA(I,IP)*X(J)*S(N7) LSQLS2.445
630 R(I)=B(I)-SJ LSQLS2.446
IF (ITMAX.LE.0) GO TO 750 LSQLS2.447
C LSQLS2.448
C CALC. SUP NORM OF X. LSQLS2.449
C LSQLS2.450
TOP=0.E0 LSQLS2.451
DO 640 J=1,N LSQLS2.452
N2=J7+J LSQLS2.453
640 TOP=MAX(TOP,ABS(X(J))*S(N2)) LSQLS2.454
C LSQLS2.455
C COMPARE RELATIVE CHANGE IN X WITH TOLERANCE EPS . LSQLS2.456
C LSQLS2.457
IF (TOP1-TOP*EPS2) 690,650,650 LSQLS2.458
650 IF (IT-ITMAX) 660,680,680 LSQLS2.459
660 IT=IT+1 LSQLS2.460
IF (IT.EQ.1) GO TO 670 LSQLS2.461
IF (TOP1.GT..25*TOP2) GO TO 690 LSQLS2.462
670 TOP2=TOP1 LSQLS2.463
GO TO (390,460), ISW LSQLS2.464
680 IT=0 LSQLS2.465
*IF -DEF,T3E OSI1F402.73
690 SJ=0.D0 LSQLS2.466
*ELSE OSI1F402.74
690 SJ=0.0 OSI1F402.75
*ENDIF OSI1F402.76
DO 700 J=1,K LSQLS2.467
SJ=SJ+R(J)**2 LSQLS2.468
700 CONTINUE LSQLS2.469
ENORM=SQRT(SJ) LSQLS2.470
IF (IRANK.EQ.N.AND.ISW.EQ.1) GO TO 710 LSQLS2.471
GO TO 730 LSQLS2.472
710 ENM1=ENORM LSQLS2.473
C LSQLS2.474
C SAVE X ARRAY. LSQLS2.475
C LSQLS2.476
DO 720 J=1,N LSQLS2.477
N1=K+J LSQLS2.478
720 R(N1)=X(J) LSQLS2.479
ISW=2 LSQLS2.480
IT=0 LSQLS2.481
GO TO 460 LSQLS2.482
C LSQLS2.483
C CHOOSE BEST SOLUTION LSQLS2.484
C LSQLS2.485
730 IF (IRANK.LT.N) GO TO 750 LSQLS2.486
IF (ENORM.LE.ENM1) GO TO 750 LSQLS2.487
DO 740 J=1,N LSQLS2.488
N1=K+J LSQLS2.489
740 X(J)=R(N1) LSQLS2.490
ENORM=ENM1 LSQLS2.491
C LSQLS2.492
C NORM OF AX - B LOCATED IN THE CELL ENORM . LSQLS2.493
C LSQLS2.494
C LSQLS2.495
C REARRANGE VARIABLES. LSQLS2.496
C LSQLS2.497
750 DO 760 J=1,N LSQLS2.498
N1=J1+J LSQLS2.499
760 S(J)=S(N1) LSQLS2.500
DO 790 J=1,N LSQLS2.501
DO 770 I=J,N LSQLS2.502
IP=S(I) LSQLS2.503
IF (J.EQ.IP) GO TO 780 LSQLS2.504
770 CONTINUE LSQLS2.505
780 S(I)=S(J) LSQLS2.506
S(J)=J LSQLS2.507
SJ=X(J) LSQLS2.508
X(J)=X(I) LSQLS2.509
790 X(I)=SJ LSQLS2.510
C LSQLS2.511
C SCALE VARIABLES. LSQLS2.512
C LSQLS2.513
DO 800 J=1,N LSQLS2.514
N2=J2+J LSQLS2.515
800 X(J)=X(J)*S(N2) LSQLS2.516
RETURN LSQLS2.517
C LSQLS2.518
C RESTORE A. LSQLS2.519
C LSQLS2.520
810 DO 820 J=1,N LSQLS2.521
N2=J2+J LSQLS2.522
DO 820 I=1,K LSQLS2.523
820 A(I,J)=AA(I,J) LSQLS2.524
RETURN LSQLS2.525
C LSQLS2.526
C GENERATE SOLUTIONS TO THE HOMOGENEOUS EQUATION AX = 0. LSQLS2.527
C LSQLS2.528
830 IF (IRANK.EQ.N) RETURN LSQLS2.529
NS=N-IRANK LSQLS2.530
DO 840 I=1,N LSQLS2.531
DO 840 J=1,NS LSQLS2.532
840 H(I,J)=0.E0 LSQLS2.533
DO 850 J=1,NS LSQLS2.534
N2=IRANK+J LSQLS2.535
850 H(N2,J)=1.E0 LSQLS2.536
IF (IRANK.EQ.0) RETURN LSQLS2.537
DO 870 J=1,IRANK LSQLS2.538
DO 870 I=1,NS LSQLS2.539
N7=K+N+J LSQLS2.540
SJ=R(N7)*H(J,I) LSQLS2.541
DO 860 K1=IRP1,N LSQLS2.542
860 SJ=SJ+H(K1,I)*A(J,K1) LSQLS2.543
N6=K+J LSQLS2.544
BP=R(N6) LSQLS2.545
DP=BP*R(N7)*SJ LSQLS2.546
A1=DP LSQLS2.547
A2=DP-A1 LSQLS2.548
H(J,I)=H(J,I)-(A1+2.*A2) LSQLS2.549
DO 870 K1=IRP1,N LSQLS2.550
DP=BP*A(J,K1)*SJ LSQLS2.551
A1=DP LSQLS2.552
A2=DP-A1 LSQLS2.553
870 H(K1,I)=H(K1,I)-(A1+2.*A2) LSQLS2.554
C LSQLS2.555
C REARRANGE ROWS OF SOLUTION MATRIX. LSQLS2.556
C LSQLS2.557
DO 880 J=1,N LSQLS2.558
N1=J1+J LSQLS2.559
880 S(J)=S(N1) LSQLS2.560
DO 910 J=1,N LSQLS2.561
DO 890 I=J,N LSQLS2.562
IP=S(I) LSQLS2.563
IF (J.EQ.IP) GO TO 900 LSQLS2.564
890 CONTINUE LSQLS2.565
900 S(I)=S(J) LSQLS2.566
S(J)=J LSQLS2.567
DO 910 K1=1,NS LSQLS2.568
A1=H(J,K1) LSQLS2.569
H(J,K1)=H(I,K1) LSQLS2.570
910 H(I,K1)=A1 LSQLS2.571
RETURN LSQLS2.572
C LSQLS2.573
1140 FORMAT (31H0WARNING. IRANK HAS BEEN SET TO,I4,6H BUT(,1PE10.3,9H) LSQLS2.574
1 RANK IS,I4,25H. IRANK IS NOW TAKEN AS ,I4,1H.) LSQLS2.575
END LSQLS2.576
*ENDIF @DYALLOC.4437