← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/ssvdc.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/ssvdc.f
 
1
*DECK SSVDC
 
2
      SUBROUTINE SSVDC (X, LDX, N, P, S, E, U, LDU, V, LDV, WORK, JOB,
 
3
     +   INFO)
 
4
C***BEGIN PROLOGUE  SSVDC
 
5
C***PURPOSE  Perform the singular value decomposition of a rectangular
 
6
C            matrix.
 
7
C***LIBRARY   SLATEC (LINPACK)
 
8
C***CATEGORY  D6
 
9
C***TYPE      SINGLE PRECISION (SSVDC-S, DSVDC-D, CSVDC-C)
 
10
C***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX,
 
11
C             SINGULAR VALUE DECOMPOSITION
 
12
C***AUTHOR  Stewart, G. W., (U. of Maryland)
 
13
C***DESCRIPTION
 
14
C
 
15
C     SSVDC is a subroutine to reduce a real NxP matrix X by orthogonal
 
16
C     transformations U and V to diagonal form.  The elements S(I) are
 
17
C     the singular values of X.  The columns of U are the corresponding
 
18
C     left singular vectors, and the columns of V the right singular
 
19
C     vectors.
 
20
C
 
21
C     On Entry
 
22
C
 
23
C         X         REAL(LDX,P), where LDX .GE. N.
 
24
C                   X contains the matrix whose singular value
 
25
C                   decomposition is to be computed.  X is
 
26
C                   destroyed by SSVDC.
 
27
C
 
28
C         LDX       INTEGER
 
29
C                   LDX is the leading dimension of the array X.
 
30
C
 
31
C         N         INTEGER
 
32
C                   N is the number of rows of the matrix X.
 
33
C
 
34
C         P         INTEGER
 
35
C                   P is the number of columns of the matrix X.
 
36
C
 
37
C         LDU       INTEGER
 
38
C                   LDU is the leading dimension of the array U.
 
39
C                   (See below).
 
40
C
 
41
C         LDV       INTEGER
 
42
C                   LDV is the leading dimension of the array V.
 
43
C                   (See below).
 
44
C
 
45
C         WORK      REAL(N)
 
46
C                   work is a scratch array.
 
47
C
 
48
C         JOB       INTEGER
 
49
C                   JOB controls the computation of the singular
 
50
C                   vectors.  It has the decimal expansion AB
 
51
C                   with the following meaning
 
52
C
 
53
C                        A .EQ. 0  Do not compute the left singular
 
54
C                                  vectors.
 
55
C                        A .EQ. 1  Return the N left singular vectors
 
56
C                                  in U.
 
57
C                        A .GE. 2  Return the first MIN(N,P) singular
 
58
C                                  vectors in U.
 
59
C                        B .EQ. 0  Do not compute the right singular
 
60
C                                  vectors.
 
61
C                        B .EQ. 1  Return the right singular vectors
 
62
C                                  in V.
 
63
C
 
64
C     On Return
 
65
C
 
66
C         S         REAL(MM), where MM=MIN(N+1,P).
 
67
C                   The first MIN(N,P) entries of S contain the
 
68
C                   singular values of X arranged in descending
 
69
C                   order of magnitude.
 
70
C
 
71
C         E         REAL(P).
 
72
C                   E ordinarily contains zeros.  However, see the
 
73
C                   discussion of INFO for exceptions.
 
74
C
 
75
C         U         REAL(LDU,K), where LDU .GE. N.  If JOBA .EQ. 1, then
 
76
C                                   K .EQ. N.  If JOBA .GE. 2 , then
 
77
C                                   K .EQ. MIN(N,P).
 
78
C                   U contains the matrix of right singular vectors.
 
79
C                   U is not referenced if JOBA .EQ. 0.  If N .LE. P
 
80
C                   or if JOBA .EQ. 2, then U may be identified with X
 
81
C                   in the subroutine call.
 
82
C
 
83
C         V         REAL(LDV,P), where LDV .GE. P.
 
84
C                   V contains the matrix of right singular vectors.
 
85
C                   V is not referenced if JOB .EQ. 0.  If P .LE. N,
 
86
C                   then V may be identified with X in the
 
87
C                   subroutine call.
 
88
C
 
89
C         INFO      INTEGER.
 
90
C                   the singular values (and their corresponding
 
91
C                   singular vectors) S(INFO+1),S(INFO+2),...,S(M)
 
92
C                   are correct (here M=MIN(N,P)).  Thus if
 
93
C                   INFO .EQ. 0, all the singular values and their
 
94
C                   vectors are correct.  In any event, the matrix
 
95
C                   B = TRANS(U)*X*V is the bidiagonal matrix
 
96
C                   with the elements of S on its diagonal and the
 
97
C                   elements of E on its super-diagonal (TRANS(U)
 
98
C                   is the transpose of U).  Thus the singular
 
99
C                   values of X and B are the same.
 
100
C
 
101
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
 
102
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
 
103
C***ROUTINES CALLED  SAXPY, SDOT, SNRM2, SROT, SROTG, SSCAL, SSWAP
 
104
C***REVISION HISTORY  (YYMMDD)
 
105
C   790319  DATE WRITTEN
 
106
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
107
C   890531  REVISION DATE from Version 3.2
 
108
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
109
C   900326  Removed duplicate information from DESCRIPTION section.
 
110
C           (WRB)
 
111
C   920501  Reformatted the REFERENCES section.  (WRB)
 
112
C***END PROLOGUE  SSVDC
 
113
      INTEGER LDX,N,P,LDU,LDV,JOB,INFO
 
114
      REAL X(LDX,*),S(*),E(*),U(LDU,*),V(LDV,*),WORK(*)
 
115
C
 
116
C
 
117
      INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
 
118
     1        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
 
119
      REAL SDOT,T
 
120
      REAL B,C,CS,EL,EMM1,F,G,SNRM2,SCALE,SHIFT,SL,SM,SN,SMM1,T1,TEST,
 
121
     1     ZTEST
 
122
      LOGICAL WANTU,WANTV
 
123
C***FIRST EXECUTABLE STATEMENT  SSVDC
 
124
C
 
125
C     SET THE MAXIMUM NUMBER OF ITERATIONS.
 
126
C
 
127
      MAXIT = 30
 
128
C
 
129
C     DETERMINE WHAT IS TO BE COMPUTED.
 
130
C
 
131
      WANTU = .FALSE.
 
132
      WANTV = .FALSE.
 
133
      JOBU = MOD(JOB,100)/10
 
134
      NCU = N
 
135
      IF (JOBU .GT. 1) NCU = MIN(N,P)
 
136
      IF (JOBU .NE. 0) WANTU = .TRUE.
 
137
      IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
 
138
C
 
139
C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
 
140
C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
 
141
C
 
142
      INFO = 0
 
143
      NCT = MIN(N-1,P)
 
144
      NRT = MAX(0,MIN(P-2,N))
 
145
      LU = MAX(NCT,NRT)
 
146
      IF (LU .LT. 1) GO TO 170
 
147
      DO 160 L = 1, LU
 
148
         LP1 = L + 1
 
149
         IF (L .GT. NCT) GO TO 20
 
150
C
 
151
C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
 
152
C           PLACE THE L-TH DIAGONAL IN S(L).
 
153
C
 
154
            S(L) = SNRM2(N-L+1,X(L,L),1)
 
155
            IF (S(L) .EQ. 0.0E0) GO TO 10
 
156
               IF (X(L,L) .NE. 0.0E0) S(L) = SIGN(S(L),X(L,L))
 
157
               CALL SSCAL(N-L+1,1.0E0/S(L),X(L,L),1)
 
158
               X(L,L) = 1.0E0 + X(L,L)
 
159
   10       CONTINUE
 
160
            S(L) = -S(L)
 
161
   20    CONTINUE
 
162
         IF (P .LT. LP1) GO TO 50
 
163
         DO 40 J = LP1, P
 
164
            IF (L .GT. NCT) GO TO 30
 
165
            IF (S(L) .EQ. 0.0E0) GO TO 30
 
166
C
 
167
C              APPLY THE TRANSFORMATION.
 
168
C
 
169
               T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
 
170
               CALL SAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
 
171
   30       CONTINUE
 
172
C
 
173
C           PLACE THE L-TH ROW OF X INTO  E FOR THE
 
174
C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
 
175
C
 
176
            E(J) = X(L,J)
 
177
   40    CONTINUE
 
178
   50    CONTINUE
 
179
         IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
 
180
C
 
181
C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
 
182
C           MULTIPLICATION.
 
183
C
 
184
            DO 60 I = L, N
 
185
               U(I,L) = X(I,L)
 
186
   60       CONTINUE
 
187
   70    CONTINUE
 
188
         IF (L .GT. NRT) GO TO 150
 
189
C
 
190
C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
 
191
C           L-TH SUPER-DIAGONAL IN E(L).
 
192
C
 
193
            E(L) = SNRM2(P-L,E(LP1),1)
 
194
            IF (E(L) .EQ. 0.0E0) GO TO 80
 
195
               IF (E(LP1) .NE. 0.0E0) E(L) = SIGN(E(L),E(LP1))
 
196
               CALL SSCAL(P-L,1.0E0/E(L),E(LP1),1)
 
197
               E(LP1) = 1.0E0 + E(LP1)
 
198
   80       CONTINUE
 
199
            E(L) = -E(L)
 
200
            IF (LP1 .GT. N .OR. E(L) .EQ. 0.0E0) GO TO 120
 
201
C
 
202
C              APPLY THE TRANSFORMATION.
 
203
C
 
204
               DO 90 I = LP1, N
 
205
                  WORK(I) = 0.0E0
 
206
   90          CONTINUE
 
207
               DO 100 J = LP1, P
 
208
                  CALL SAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
 
209
  100          CONTINUE
 
210
               DO 110 J = LP1, P
 
211
                  CALL SAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
 
212
  110          CONTINUE
 
213
  120       CONTINUE
 
214
            IF (.NOT.WANTV) GO TO 140
 
215
C
 
216
C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
 
217
C              BACK MULTIPLICATION.
 
218
C
 
219
               DO 130 I = LP1, P
 
220
                  V(I,L) = E(I)
 
221
  130          CONTINUE
 
222
  140       CONTINUE
 
223
  150    CONTINUE
 
224
  160 CONTINUE
 
225
  170 CONTINUE
 
226
C
 
227
C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
 
228
C
 
229
      M = MIN(P,N+1)
 
230
      NCTP1 = NCT + 1
 
231
      NRTP1 = NRT + 1
 
232
      IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
 
233
      IF (N .LT. M) S(M) = 0.0E0
 
234
      IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
 
235
      E(M) = 0.0E0
 
236
C
 
237
C     IF REQUIRED, GENERATE U.
 
238
C
 
239
      IF (.NOT.WANTU) GO TO 300
 
240
         IF (NCU .LT. NCTP1) GO TO 200
 
241
         DO 190 J = NCTP1, NCU
 
242
            DO 180 I = 1, N
 
243
               U(I,J) = 0.0E0
 
244
  180       CONTINUE
 
245
            U(J,J) = 1.0E0
 
246
  190    CONTINUE
 
247
  200    CONTINUE
 
248
         IF (NCT .LT. 1) GO TO 290
 
249
         DO 280 LL = 1, NCT
 
250
            L = NCT - LL + 1
 
251
            IF (S(L) .EQ. 0.0E0) GO TO 250
 
252
               LP1 = L + 1
 
253
               IF (NCU .LT. LP1) GO TO 220
 
254
               DO 210 J = LP1, NCU
 
255
                  T = -SDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
 
256
                  CALL SAXPY(N-L+1,T,U(L,L),1,U(L,J),1)
 
257
  210          CONTINUE
 
258
  220          CONTINUE
 
259
               CALL SSCAL(N-L+1,-1.0E0,U(L,L),1)
 
260
               U(L,L) = 1.0E0 + U(L,L)
 
261
               LM1 = L - 1
 
262
               IF (LM1 .LT. 1) GO TO 240
 
263
               DO 230 I = 1, LM1
 
264
                  U(I,L) = 0.0E0
 
265
  230          CONTINUE
 
266
  240          CONTINUE
 
267
            GO TO 270
 
268
  250       CONTINUE
 
269
               DO 260 I = 1, N
 
270
                  U(I,L) = 0.0E0
 
271
  260          CONTINUE
 
272
               U(L,L) = 1.0E0
 
273
  270       CONTINUE
 
274
  280    CONTINUE
 
275
  290    CONTINUE
 
276
  300 CONTINUE
 
277
C
 
278
C     IF IT IS REQUIRED, GENERATE V.
 
279
C
 
280
      IF (.NOT.WANTV) GO TO 350
 
281
         DO 340 LL = 1, P
 
282
            L = P - LL + 1
 
283
            LP1 = L + 1
 
284
            IF (L .GT. NRT) GO TO 320
 
285
            IF (E(L) .EQ. 0.0E0) GO TO 320
 
286
               DO 310 J = LP1, P
 
287
                  T = -SDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
 
288
                  CALL SAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1)
 
289
  310          CONTINUE
 
290
  320       CONTINUE
 
291
            DO 330 I = 1, P
 
292
               V(I,L) = 0.0E0
 
293
  330       CONTINUE
 
294
            V(L,L) = 1.0E0
 
295
  340    CONTINUE
 
296
  350 CONTINUE
 
297
C
 
298
C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
 
299
C
 
300
      MM = M
 
301
      ITER = 0
 
302
  360 CONTINUE
 
303
C
 
304
C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
 
305
C
 
306
         IF (M .EQ. 0) GO TO 620
 
307
C
 
308
C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
 
309
C        FLAG AND RETURN.
 
310
C
 
311
         IF (ITER .LT. MAXIT) GO TO 370
 
312
            INFO = M
 
313
            GO TO 620
 
314
  370    CONTINUE
 
315
C
 
316
C        THIS SECTION OF THE PROGRAM INSPECTS FOR
 
317
C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
 
318
C        COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
 
319
C
 
320
C           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
 
321
C           KASE = 2     IF S(L) IS NEGLIGIBLE AND L.LT.M
 
322
C           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
 
323
C                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
 
324
C           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
 
325
C
 
326
         DO 390 LL = 1, M
 
327
            L = M - LL
 
328
            IF (L .EQ. 0) GO TO 400
 
329
            TEST = ABS(S(L)) + ABS(S(L+1))
 
330
            ZTEST = TEST + ABS(E(L))
 
331
            IF (ZTEST .NE. TEST) GO TO 380
 
332
               E(L) = 0.0E0
 
333
               GO TO 400
 
334
  380       CONTINUE
 
335
  390    CONTINUE
 
336
  400    CONTINUE
 
337
         IF (L .NE. M - 1) GO TO 410
 
338
            KASE = 4
 
339
         GO TO 480
 
340
  410    CONTINUE
 
341
            LP1 = L + 1
 
342
            MP1 = M + 1
 
343
            DO 430 LLS = LP1, MP1
 
344
               LS = M - LLS + LP1
 
345
               IF (LS .EQ. L) GO TO 440
 
346
               TEST = 0.0E0
 
347
               IF (LS .NE. M) TEST = TEST + ABS(E(LS))
 
348
               IF (LS .NE. L + 1) TEST = TEST + ABS(E(LS-1))
 
349
               ZTEST = TEST + ABS(S(LS))
 
350
               IF (ZTEST .NE. TEST) GO TO 420
 
351
                  S(LS) = 0.0E0
 
352
                  GO TO 440
 
353
  420          CONTINUE
 
354
  430       CONTINUE
 
355
  440       CONTINUE
 
356
            IF (LS .NE. L) GO TO 450
 
357
               KASE = 3
 
358
            GO TO 470
 
359
  450       CONTINUE
 
360
            IF (LS .NE. M) GO TO 460
 
361
               KASE = 1
 
362
            GO TO 470
 
363
  460       CONTINUE
 
364
               KASE = 2
 
365
               L = LS
 
366
  470       CONTINUE
 
367
  480    CONTINUE
 
368
         L = L + 1
 
369
C
 
370
C        PERFORM THE TASK INDICATED BY KASE.
 
371
C
 
372
         GO TO (490,520,540,570), KASE
 
373
C
 
374
C        DEFLATE NEGLIGIBLE S(M).
 
375
C
 
376
  490    CONTINUE
 
377
            MM1 = M - 1
 
378
            F = E(M-1)
 
379
            E(M-1) = 0.0E0
 
380
            DO 510 KK = L, MM1
 
381
               K = MM1 - KK + L
 
382
               T1 = S(K)
 
383
               CALL SROTG(T1,F,CS,SN)
 
384
               S(K) = T1
 
385
               IF (K .EQ. L) GO TO 500
 
386
                  F = -SN*E(K-1)
 
387
                  E(K-1) = CS*E(K-1)
 
388
  500          CONTINUE
 
389
               IF (WANTV) CALL SROT(P,V(1,K),1,V(1,M),1,CS,SN)
 
390
  510       CONTINUE
 
391
         GO TO 610
 
392
C
 
393
C        SPLIT AT NEGLIGIBLE S(L).
 
394
C
 
395
  520    CONTINUE
 
396
            F = E(L-1)
 
397
            E(L-1) = 0.0E0
 
398
            DO 530 K = L, M
 
399
               T1 = S(K)
 
400
               CALL SROTG(T1,F,CS,SN)
 
401
               S(K) = T1
 
402
               F = -SN*E(K)
 
403
               E(K) = CS*E(K)
 
404
               IF (WANTU) CALL SROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
 
405
  530       CONTINUE
 
406
         GO TO 610
 
407
C
 
408
C        PERFORM ONE QR STEP.
 
409
C
 
410
  540    CONTINUE
 
411
C
 
412
C           CALCULATE THE SHIFT.
 
413
C
 
414
            SCALE = MAX(ABS(S(M)),ABS(S(M-1)),ABS(E(M-1)),ABS(S(L)),
 
415
     1                    ABS(E(L)))
 
416
            SM = S(M)/SCALE
 
417
            SMM1 = S(M-1)/SCALE
 
418
            EMM1 = E(M-1)/SCALE
 
419
            SL = S(L)/SCALE
 
420
            EL = E(L)/SCALE
 
421
            B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0E0
 
422
            C = (SM*EMM1)**2
 
423
            SHIFT = 0.0E0
 
424
            IF (B .EQ. 0.0E0 .AND. C .EQ. 0.0E0) GO TO 550
 
425
               SHIFT = SQRT(B**2+C)
 
426
               IF (B .LT. 0.0E0) SHIFT = -SHIFT
 
427
               SHIFT = C/(B + SHIFT)
 
428
  550       CONTINUE
 
429
            F = (SL + SM)*(SL - SM) - SHIFT
 
430
            G = SL*EL
 
431
C
 
432
C           CHASE ZEROS.
 
433
C
 
434
            MM1 = M - 1
 
435
            DO 560 K = L, MM1
 
436
               CALL SROTG(F,G,CS,SN)
 
437
               IF (K .NE. L) E(K-1) = F
 
438
               F = CS*S(K) + SN*E(K)
 
439
               E(K) = CS*E(K) - SN*S(K)
 
440
               G = SN*S(K+1)
 
441
               S(K+1) = CS*S(K+1)
 
442
               IF (WANTV) CALL SROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
 
443
               CALL SROTG(F,G,CS,SN)
 
444
               S(K) = F
 
445
               F = CS*E(K) + SN*S(K+1)
 
446
               S(K+1) = -SN*E(K) + CS*S(K+1)
 
447
               G = SN*E(K+1)
 
448
               E(K+1) = CS*E(K+1)
 
449
               IF (WANTU .AND. K .LT. N)
 
450
     1            CALL SROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
 
451
  560       CONTINUE
 
452
            E(M-1) = F
 
453
            ITER = ITER + 1
 
454
         GO TO 610
 
455
C
 
456
C        CONVERGENCE.
 
457
C
 
458
  570    CONTINUE
 
459
C
 
460
C           MAKE THE SINGULAR VALUE  POSITIVE.
 
461
C
 
462
            IF (S(L) .GE. 0.0E0) GO TO 580
 
463
               S(L) = -S(L)
 
464
               IF (WANTV) CALL SSCAL(P,-1.0E0,V(1,L),1)
 
465
  580       CONTINUE
 
466
C
 
467
C           ORDER THE SINGULAR VALUE.
 
468
C
 
469
  590       IF (L .EQ. MM) GO TO 600
 
470
               IF (S(L) .GE. S(L+1)) GO TO 600
 
471
               T = S(L)
 
472
               S(L) = S(L+1)
 
473
               S(L+1) = T
 
474
               IF (WANTV .AND. L .LT. P)
 
475
     1            CALL SSWAP(P,V(1,L),1,V(1,L+1),1)
 
476
               IF (WANTU .AND. L .LT. N)
 
477
     1            CALL SSWAP(N,U(1,L),1,U(1,L+1),1)
 
478
               L = L + 1
 
479
            GO TO 590
 
480
  600       CONTINUE
 
481
            ITER = 0
 
482
            M = M - 1
 
483
  610    CONTINUE
 
484
      GO TO 360
 
485
  620 CONTINUE
 
486
      RETURN
 
487
      END