← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/la05as.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/la05as.f
 
1
*DECK LA05AS
 
2
      SUBROUTINE LA05AS (A, IND, NZ, IA, N, IP, IW, W, G, U)
 
3
C***BEGIN PROLOGUE  LA05AS
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Subsidiary to SPLP
 
6
C***LIBRARY   SLATEC
 
7
C***TYPE      SINGLE PRECISION (LA05AS-S, LA05AD-D)
 
8
C***AUTHOR  (UNKNOWN)
 
9
C***DESCRIPTION
 
10
C
 
11
C     THIS SUBPROGRAM IS A SLIGHT MODIFICATION OF A SUBPROGRAM
 
12
C     FROM THE C. 1979 AERE HARWELL LIBRARY.  THE NAME OF THE
 
13
C     CORRESPONDING HARWELL CODE CAN BE OBTAINED BY DELETING
 
14
C     THE FINAL LETTER =S= IN THE NAMES USED HERE.
 
15
C     REVISIONS MADE BY R J HANSON, SNLA, AUGUST, 1979.
 
16
C     REVISED SEP. 13, 1979.
 
17
C
 
18
C     ROYALTIES HAVE BEEN PAID TO AERE-UK FOR USE OF THEIR CODES
 
19
C     IN THE PACKAGE GIVEN HERE.  ANY PRIMARY USAGE OF THE HARWELL
 
20
C     SUBROUTINES REQUIRES A ROYALTY AGREEMENT AND PAYMENT BETWEEN
 
21
C     THE USER AND AERE-UK.  ANY USAGE OF THE SANDIA WRITTEN CODES
 
22
C     SPLP( ) (WHICH USES THE HARWELL SUBROUTINES) IS PERMITTED.
 
23
C
 
24
C IP(I,1),IP(I,2) POINT TO THE START OF ROW/COL I.
 
25
C IW(I,1),IW(I,2) HOLD THE NUMBER OF NON-ZEROS IN ROW/COL I.
 
26
C DURING THE MAIN BODY OF THIS SUBROUTINE THE VECTORS IW(.,3),IW(.,5),
 
27
C     IW(.,7) ARE USED TO HOLD DOUBLY LINKED LISTS OF ROWS THAT HAVE
 
28
C     NOT BEEN PIVOTAL AND HAVE EQUAL NUMBERS OF NON-ZEROS.
 
29
C IW(.,4),IW(.,6),IW(.,8) HOLD SIMILAR LISTS FOR THE COLUMNS.
 
30
C IW(I,3),IW(I,4) HOLD FIRST ROW/COLUMN TO HAVE I NON-ZEROS
 
31
C     OR ZERO IF THERE ARE NONE.
 
32
C IW(I,5), IW(I,6) HOLD ROW/COL NUMBER OF ROW/COL PRIOR TO ROW/COL I
 
33
C     IN ITS LIST, OR ZERO IF NONE.
 
34
C IW(I,7), IW(I,8) HOLD ROW/COL NUMBER OF ROW/COL AFTER ROW/COL I
 
35
C     IN ITS LIST, OR ZERO IF NONE.
 
36
C FOR ROWS/COLS THAT HAVE BEEN PIVOTAL IW(I,5),IW(I,6) HOLD NEGATION OF
 
37
C     POSITION OF ROW/COL I IN THE PIVOTAL ORDERING.
 
38
C
 
39
C***SEE ALSO  SPLP
 
40
C***ROUTINES CALLED  LA05ES, MC20AS, R1MACH, XERMSG, XSETUN
 
41
C***COMMON BLOCKS    LA05DS
 
42
C***REVISION HISTORY  (YYMMDD)
 
43
C   811215  DATE WRITTEN
 
44
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
45
C   890605  Corrected references to XERRWV.  (WRB)
 
46
C   890831  Modified array declarations.  (WRB)
 
47
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
48
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 
49
C   900402  Added TYPE section.  (WRB)
 
50
C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
 
51
C***END PROLOGUE  LA05AS
 
52
      INTEGER IP(N,2)
 
53
      INTEGER IND(IA,2), IW(N,8)
 
54
      REAL A(*), AMAX, AU, AM, G, U, SMALL, W(*)
 
55
      LOGICAL FIRST
 
56
      CHARACTER*8 XERN0, XERN1, XERN2
 
57
C
 
58
      COMMON /LA05DS/ SMALL, LP, LENL, LENU, NCP, LROW, LCOL
 
59
C EPS IS THE RELATIVE ACCURACY OF FLOATING-POINT COMPUTATION
 
60
      SAVE EPS, FIRST
 
61
      DATA FIRST /.TRUE./
 
62
C***FIRST EXECUTABLE STATEMENT  LA05AS
 
63
      IF (FIRST) THEN
 
64
         EPS = 2.0E0 * R1MACH(4)
 
65
      ENDIF
 
66
      FIRST = .FALSE.
 
67
C
 
68
C     SET THE OUTPUT UNIT NUMBER FOR THE ERROR PROCESSOR.
 
69
C     THE USAGE OF THIS ERROR PROCESSOR IS DOCUMENTED IN THE
 
70
C     SANDIA LABS. TECH. REPT. SAND78-1189, BY R E JONES.
 
71
      CALL XSETUN(LP)
 
72
      IF (U.GT.1.0E0) U = 1.0E0
 
73
      IF (U.LT.EPS) U = EPS
 
74
      IF (N.LT.1) GO TO 670
 
75
      G = 0.
 
76
      DO 50 I=1,N
 
77
         W(I) = 0.
 
78
         DO 40 J=1,5
 
79
            IW(I,J) = 0
 
80
   40    CONTINUE
 
81
   50 CONTINUE
 
82
C
 
83
C FLUSH OUT SMALL ENTRIES, COUNT ELEMENTS IN ROWS AND COLUMNS
 
84
      L = 1
 
85
      LENU = NZ
 
86
      DO 80 IDUMMY=1,NZ
 
87
         IF (L.GT.LENU) GO TO 90
 
88
         DO 60 K=L,LENU
 
89
            IF (ABS(A(K)).LE.SMALL) GO TO 70
 
90
            I = IND(K,1)
 
91
            J = IND(K,2)
 
92
            G = MAX(ABS(A(K)),G)
 
93
            IF (I.LT.1 .OR. I.GT.N) GO TO 680
 
94
            IF (J.LT.1 .OR. J.GT.N) GO TO 680
 
95
            IW(I,1) = IW(I,1) + 1
 
96
            IW(J,2) = IW(J,2) + 1
 
97
   60    CONTINUE
 
98
         GO TO 90
 
99
   70    L = K
 
100
         A(L) = A(LENU)
 
101
         IND(L,1) = IND(LENU,1)
 
102
         IND(L,2) = IND(LENU,2)
 
103
         LENU = LENU - 1
 
104
   80 CONTINUE
 
105
C
 
106
   90 LENL = 0
 
107
      LROW = LENU
 
108
      LCOL = LROW
 
109
C MCP IS THE MAXIMUM NUMBER OF COMPRESSES PERMITTED BEFORE AN
 
110
C     ERROR RETURN RESULTS.
 
111
      MCP = MAX(N/10,20)
 
112
      NCP = 0
 
113
C CHECK FOR NULL ROW OR COLUMN AND INITIALIZE IP(I,2) TO POINT
 
114
C     JUST BEYOND WHERE THE LAST COMPONENT OF COLUMN I OF A WILL
 
115
C     BE STORED.
 
116
      K = 1
 
117
      DO 110 IR=1,N
 
118
         K = K + IW(IR,2)
 
119
         IP(IR,2) = K
 
120
         DO 100 L=1,2
 
121
            IF (IW(IR,L).LE.0) GO TO 700
 
122
  100    CONTINUE
 
123
  110 CONTINUE
 
124
C REORDER BY ROWS
 
125
C CHECK FOR DOUBLE ENTRIES WHILE USING THE NEWLY CONSTRUCTED
 
126
C     ROW FILE TO CONSTRUCT THE COLUMN FILE. NOTE THAT BY PUTTING
 
127
C    THE ENTRIES IN BACKWARDS AND DECREASING IP(J,2) EACH TIME IT
 
128
C     IS USED WE AUTOMATICALLY LEAVE IT POINTING TO THE FIRST ELEMENT.
 
129
      CALL MC20AS(N, LENU, A, IND(1,2), IP, IND(1,1), 0)
 
130
      KL = LENU
 
131
      DO 130 II=1,N
 
132
         IR = N + 1 - II
 
133
         KP = IP(IR,1)
 
134
         DO 120 K=KP,KL
 
135
            J = IND(K,2)
 
136
            IF (IW(J,5).EQ.IR) GO TO 660
 
137
            IW(J,5) = IR
 
138
            KR = IP(J,2) - 1
 
139
            IP(J,2) = KR
 
140
            IND(KR,1) = IR
 
141
  120    CONTINUE
 
142
         KL = KP - 1
 
143
  130 CONTINUE
 
144
C
 
145
C SET UP LINKED LISTS OF ROWS AND COLS WITH EQUAL NUMBERS OF NON-ZEROS.
 
146
      DO 150 L=1,2
 
147
         DO 140 I=1,N
 
148
            NZ = IW(I,L)
 
149
            IN = IW(NZ,L+2)
 
150
            IW(NZ,L+2) = I
 
151
            IW(I,L+6) = IN
 
152
            IW(I,L+4) = 0
 
153
            IF (IN.NE.0) IW(IN,L+4) = I
 
154
  140    CONTINUE
 
155
  150 CONTINUE
 
156
C
 
157
C
 
158
C START OF MAIN ELIMINATION LOOP.
 
159
      DO 590 IPV=1,N
 
160
C FIND PIVOT. JCOST IS MARKOWITZ COST OF CHEAPEST PIVOT FOUND SO FAR,
 
161
C     WHICH IS IN ROW IPP AND COLUMN JP.
 
162
         JCOST = N*N
 
163
C LOOP ON LENGTH OF COLUMN TO BE SEARCHED
 
164
         DO 240 NZ=1,N
 
165
            IF (JCOST.LE.(NZ-1)**2) GO TO 250
 
166
            J = IW(NZ,4)
 
167
C SEARCH COLUMNS WITH NZ NON-ZEROS.
 
168
            DO 190 IDUMMY=1,N
 
169
               IF (J.LE.0) GO TO 200
 
170
               KP = IP(J,2)
 
171
               KL = KP + IW(J,2) - 1
 
172
               DO 180 K=KP,KL
 
173
                  I = IND(K,1)
 
174
                  KCOST = (NZ-1)*(IW(I,1)-1)
 
175
                  IF (KCOST.GE.JCOST) GO TO 180
 
176
                  IF (NZ.EQ.1) GO TO 170
 
177
C FIND LARGEST ELEMENT IN ROW OF POTENTIAL PIVOT.
 
178
                  AMAX = 0.
 
179
                  K1 = IP(I,1)
 
180
                  K2 = IW(I,1) + K1 - 1
 
181
                  DO 160 KK=K1,K2
 
182
                     AMAX = MAX(AMAX,ABS(A(KK)))
 
183
                     IF (IND(KK,2).EQ.J) KJ = KK
 
184
  160             CONTINUE
 
185
C PERFORM STABILITY TEST.
 
186
                  IF (ABS(A(KJ)).LT.AMAX*U) GO TO 180
 
187
  170             JCOST = KCOST
 
188
                  IPP = I
 
189
                  JP = J
 
190
                  IF (JCOST.LE.(NZ-1)**2) GO TO 250
 
191
  180          CONTINUE
 
192
               J = IW(J,8)
 
193
  190       CONTINUE
 
194
C SEARCH ROWS WITH NZ NON-ZEROS.
 
195
  200       I = IW(NZ,3)
 
196
            DO 230 IDUMMY=1,N
 
197
               IF (I.LE.0) GO TO 240
 
198
               AMAX = 0.
 
199
               KP = IP(I,1)
 
200
               KL = KP + IW(I,1) - 1
 
201
C FIND LARGEST ELEMENT IN THE ROW
 
202
               DO 210 K=KP,KL
 
203
                  AMAX = MAX(ABS(A(K)),AMAX)
 
204
  210          CONTINUE
 
205
               AU = AMAX*U
 
206
               DO 220 K=KP,KL
 
207
C PERFORM STABILITY TEST.
 
208
                  IF (ABS(A(K)).LT.AU) GO TO 220
 
209
                  J = IND(K,2)
 
210
                  KCOST = (NZ-1)*(IW(J,2)-1)
 
211
                  IF (KCOST.GE.JCOST) GO TO 220
 
212
                  JCOST = KCOST
 
213
                  IPP = I
 
214
                  JP = J
 
215
                  IF (JCOST.LE.(NZ-1)**2) GO TO 250
 
216
  220          CONTINUE
 
217
               I = IW(I,7)
 
218
  230       CONTINUE
 
219
  240    CONTINUE
 
220
C
 
221
C PIVOT FOUND.
 
222
C REMOVE ROWS AND COLUMNS INVOLVED IN ELIMINATION FROM ORDERING VECTORS.
 
223
  250    KP = IP(JP,2)
 
224
         KL = IW(JP,2) + KP - 1
 
225
         DO 290 L=1,2
 
226
            DO 280 K=KP,KL
 
227
               I = IND(K,L)
 
228
               IL = IW(I,L+4)
 
229
               IN = IW(I,L+6)
 
230
               IF (IL.EQ.0) GO TO 260
 
231
               IW(IL,L+6) = IN
 
232
               GO TO 270
 
233
  260          NZ = IW(I,L)
 
234
               IW(NZ,L+2) = IN
 
235
  270          IF (IN.GT.0) IW(IN,L+4) = IL
 
236
  280       CONTINUE
 
237
            KP = IP(IPP,1)
 
238
            KL = KP + IW(IPP,1) - 1
 
239
  290    CONTINUE
 
240
C STORE PIVOT
 
241
         IW(IPP,5) = -IPV
 
242
         IW(JP,6) = -IPV
 
243
C ELIMINATE PIVOTAL ROW FROM COLUMN FILE AND FIND PIVOT IN ROW FILE.
 
244
         DO 320 K=KP,KL
 
245
            J = IND(K,2)
 
246
            KPC = IP(J,2)
 
247
            IW(J,2) = IW(J,2) - 1
 
248
            KLC = KPC + IW(J,2)
 
249
            DO 300 KC=KPC,KLC
 
250
               IF (IPP.EQ.IND(KC,1)) GO TO 310
 
251
  300       CONTINUE
 
252
  310       IND(KC,1) = IND(KLC,1)
 
253
            IND(KLC,1) = 0
 
254
            IF (J.EQ.JP) KR = K
 
255
  320    CONTINUE
 
256
C BRING PIVOT TO FRONT OF PIVOTAL ROW.
 
257
         AU = A(KR)
 
258
         A(KR) = A(KP)
 
259
         A(KP) = AU
 
260
         IND(KR,2) = IND(KP,2)
 
261
         IND(KP,2) = JP
 
262
C
 
263
C PERFORM ELIMINATION ITSELF, LOOPING ON NON-ZEROS IN PIVOT COLUMN.
 
264
         NZC = IW(JP,2)
 
265
         IF (NZC.EQ.0) GO TO 550
 
266
         DO 540 NC=1,NZC
 
267
            KC = IP(JP,2) + NC - 1
 
268
            IR = IND(KC,1)
 
269
C SEARCH NON-PIVOT ROW FOR ELEMENT TO BE ELIMINATED.
 
270
            KR = IP(IR,1)
 
271
            KRL = KR + IW(IR,1) - 1
 
272
            DO 330 KNP=KR,KRL
 
273
               IF (JP.EQ.IND(KNP,2)) GO TO 340
 
274
  330       CONTINUE
 
275
C BRING ELEMENT TO BE ELIMINATED TO FRONT OF ITS ROW.
 
276
  340       AM = A(KNP)
 
277
            A(KNP) = A(KR)
 
278
            A(KR) = AM
 
279
            IND(KNP,2) = IND(KR,2)
 
280
            IND(KR,2) = JP
 
281
            AM = -A(KR)/A(KP)
 
282
C COMPRESS ROW FILE UNLESS IT IS CERTAIN THAT THERE IS ROOM FOR NEW ROW.
 
283
            IF (LROW+IW(IR,1)+IW(IPP,1)+LENL.LE.IA) GO TO 350
 
284
            IF (NCP.GE.MCP .OR. LENU+IW(IR,1)+IW(IPP,1)+LENL.GT.IA) GO
 
285
     *       TO 710
 
286
            CALL LA05ES(A, IND(1,2), IP, N, IW, IA, .TRUE.)
 
287
            KP = IP(IPP,1)
 
288
            KR = IP(IR,1)
 
289
  350       KRL = KR + IW(IR,1) - 1
 
290
            KQ = KP + 1
 
291
            KPL = KP + IW(IPP,1) - 1
 
292
C PLACE PIVOT ROW (EXCLUDING PIVOT ITSELF) IN W.
 
293
            IF (KQ.GT.KPL) GO TO 370
 
294
            DO 360 K=KQ,KPL
 
295
               J = IND(K,2)
 
296
               W(J) = A(K)
 
297
  360       CONTINUE
 
298
  370       IP(IR,1) = LROW + 1
 
299
C
 
300
C TRANSFER MODIFIED ELEMENTS.
 
301
            IND(KR,2) = 0
 
302
            KR = KR + 1
 
303
            IF (KR.GT.KRL) GO TO 430
 
304
            DO 420 KS=KR,KRL
 
305
               J = IND(KS,2)
 
306
               AU = A(KS) + AM*W(J)
 
307
               IND(KS,2) = 0
 
308
C IF ELEMENT IS VERY SMALL REMOVE IT FROM U.
 
309
               IF (ABS(AU).LE.SMALL) GO TO 380
 
310
               G = MAX(G,ABS(AU))
 
311
               LROW = LROW + 1
 
312
               A(LROW) = AU
 
313
               IND(LROW,2) = J
 
314
               GO TO 410
 
315
  380          LENU = LENU - 1
 
316
C REMOVE ELEMENT FROM COL FILE.
 
317
               K = IP(J,2)
 
318
               KL = K + IW(J,2) - 1
 
319
               IW(J,2) = KL - K
 
320
               DO 390 KK=K,KL
 
321
                  IF (IND(KK,1).EQ.IR) GO TO 400
 
322
  390          CONTINUE
 
323
  400          IND(KK,1) = IND(KL,1)
 
324
               IND(KL,1) = 0
 
325
  410          W(J) = 0.
 
326
  420       CONTINUE
 
327
C
 
328
C SCAN PIVOT ROW FOR FILLS.
 
329
  430       IF (KQ.GT.KPL) GO TO 520
 
330
            DO 510 KS=KQ,KPL
 
331
               J = IND(KS,2)
 
332
               AU = AM*W(J)
 
333
               IF (ABS(AU).LE.SMALL) GO TO 500
 
334
               LROW = LROW + 1
 
335
               A(LROW) = AU
 
336
               IND(LROW,2) = J
 
337
               LENU = LENU + 1
 
338
C
 
339
C CREATE FILL IN COLUMN FILE.
 
340
               NZ = IW(J,2)
 
341
               K = IP(J,2)
 
342
               KL = K + NZ - 1
 
343
               IF (NZ .EQ. 0) GO TO 460
 
344
C IF POSSIBLE PLACE NEW ELEMENT AT END OF PRESENT ENTRY.
 
345
               IF (KL.NE.LCOL) GO TO 440
 
346
               IF (LCOL+LENL.GE.IA) GO TO 460
 
347
               LCOL = LCOL + 1
 
348
               GO TO 450
 
349
  440          IF (IND(KL+1,1).NE.0) GO TO 460
 
350
  450          IND(KL+1,1) = IR
 
351
               GO TO 490
 
352
C NEW ENTRY HAS TO BE CREATED.
 
353
  460          IF (LCOL+LENL+NZ+1.LT.IA) GO TO 470
 
354
C COMPRESS COLUMN FILE IF THERE IS NOT ROOM FOR NEW ENTRY.
 
355
               IF (NCP.GE.MCP .OR. LENU+LENL+NZ+1.GE.IA) GO TO 710
 
356
               CALL LA05ES(A, IND, IP(1,2), N, IW(1,2), IA, .FALSE.)
 
357
               K = IP(J,2)
 
358
               KL = K + NZ - 1
 
359
C TRANSFER OLD ENTRY INTO NEW.
 
360
  470          IP(J,2) = LCOL + 1
 
361
               IF (KL .LT. K) GO TO 485
 
362
               DO 480 KK=K,KL
 
363
                  LCOL = LCOL + 1
 
364
                  IND(LCOL,1) = IND(KK,1)
 
365
                  IND(KK,1) = 0
 
366
  480          CONTINUE
 
367
  485          CONTINUE
 
368
C ADD NEW ELEMENT.
 
369
               LCOL = LCOL + 1
 
370
               IND(LCOL,1) = IR
 
371
  490          G = MAX(G,ABS(AU))
 
372
               IW(J,2) = NZ + 1
 
373
  500          W(J) = 0.
 
374
  510       CONTINUE
 
375
  520       IW(IR,1) = LROW + 1 - IP(IR,1)
 
376
C
 
377
C STORE MULTIPLIER
 
378
            IF (LENL+LCOL+1.LE.IA) GO TO 530
 
379
C COMPRESS COL FILE IF NECESSARY.
 
380
            IF (NCP.GE.MCP) GO TO 710
 
381
            CALL LA05ES(A, IND, IP(1,2), N, IW(1,2), IA, .FALSE.)
 
382
  530       K = IA - LENL
 
383
            LENL = LENL + 1
 
384
            A(K) = AM
 
385
            IND(K,1) = IPP
 
386
            IND(K,2) = IR
 
387
            LENU = LENU - 1
 
388
  540    CONTINUE
 
389
C
 
390
C INSERT ROWS AND COLUMNS INVOLVED IN ELIMINATION IN LINKED LISTS
 
391
C     OF EQUAL NUMBERS OF NON-ZEROS.
 
392
  550    K1 = IP(JP,2)
 
393
         K2 = IW(JP,2) + K1 - 1
 
394
         IW(JP,2) = 0
 
395
         DO 580 L=1,2
 
396
            IF (K2.LT.K1) GO TO 570
 
397
            DO 560 K=K1,K2
 
398
               IR = IND(K,L)
 
399
               IF (L.EQ.1) IND(K,L) = 0
 
400
               NZ = IW(IR,L)
 
401
               IF (NZ.LE.0) GO TO 720
 
402
               IN = IW(NZ,L+2)
 
403
               IW(IR,L+6) = IN
 
404
               IW(IR,L+4) = 0
 
405
               IW(NZ,L+2) = IR
 
406
               IF (IN.NE.0) IW(IN,L+4) = IR
 
407
  560       CONTINUE
 
408
  570       K1 = IP(IPP,1) + 1
 
409
            K2 = IW(IPP,1) + K1 - 2
 
410
  580    CONTINUE
 
411
  590 CONTINUE
 
412
C
 
413
C RESET COLUMN FILE TO REFER TO U AND STORE ROW/COL NUMBERS IN
 
414
C     PIVOTAL ORDER IN IW(.,3),IW(.,4)
 
415
      DO 600 I=1,N
 
416
         J = -IW(I,5)
 
417
         IW(J,3) = I
 
418
         J = -IW(I,6)
 
419
         IW(J,4) = I
 
420
         IW(I,2) = 0
 
421
  600 CONTINUE
 
422
      DO 620 I=1,N
 
423
         KP = IP(I,1)
 
424
         KL = IW(I,1) + KP - 1
 
425
         DO 610 K=KP,KL
 
426
            J = IND(K,2)
 
427
            IW(J,2) = IW(J,2) + 1
 
428
  610    CONTINUE
 
429
  620 CONTINUE
 
430
      K = 1
 
431
      DO 630 I=1,N
 
432
         K = K + IW(I,2)
 
433
         IP(I,2) = K
 
434
  630 CONTINUE
 
435
      LCOL = K - 1
 
436
      DO 650 II=1,N
 
437
         I = IW(II,3)
 
438
         KP = IP(I,1)
 
439
         KL = IW(I,1) + KP - 1
 
440
         DO 640 K=KP,KL
 
441
            J = IND(K,2)
 
442
            KN = IP(J,2) - 1
 
443
            IP(J,2) = KN
 
444
            IND(KN,1) = I
 
445
  640    CONTINUE
 
446
  650 CONTINUE
 
447
      RETURN
 
448
C
 
449
C     THE FOLLOWING INSTRUCTIONS IMPLEMENT THE FAILURE EXITS.
 
450
C
 
451
  660 IF (LP.GT.0) THEN
 
452
         WRITE (XERN1, '(I8)') IR
 
453
         WRITE (XERN2, '(I8)') J
 
454
         CALL XERMSG ('SLATEC', 'LA05AS', 'MORE THAN ONE MATRIX ' //
 
455
     *      'ENTRY.  HERE ROW = ' // XERN1 // ' AND COL = ' // XERN2,
 
456
     *      -4, 1)
 
457
      ENDIF
 
458
      G = -4.
 
459
      RETURN
 
460
C
 
461
  670 IF (LP.GT.0) CALL XERMSG ('SLATEC', 'LA05AS',
 
462
     *   'THE ORDER OF THE SYSTEM, N, IS NOT POSITIVE.', -1, 1)
 
463
      G = -1.0E0
 
464
      RETURN
 
465
C
 
466
  680 IF (LP.GT.0) THEN
 
467
         WRITE (XERN0, '(I8)') K
 
468
         WRITE (XERN1, '(I8)') I
 
469
         WRITE (XERN2, '(I8)') J
 
470
         CALL XERMSG ('SLATEC', 'LA05AS', 'ELEMENT K = ' // XERN0 //
 
471
     *      ' IS OUT OF BOUNDS.$$HERE ROW = ' // XERN1 //
 
472
     *      ' AND COL = ' // XERN2, -3, 1)
 
473
      ENDIF
 
474
      G = -3.
 
475
      RETURN
 
476
C
 
477
  700 IF (LP.GT.0) THEN
 
478
         WRITE (XERN1, '(I8)') L
 
479
         CALL XERMSG ('SLATEC', 'LA05AS', 'ROW OR COLUMN HAS NO ' //
 
480
     *      'ELEMENTS.  HERE INDEX = ' // XERN1, -2, 1)
 
481
      ENDIF
 
482
      G = -2.
 
483
      RETURN
 
484
C
 
485
  710 IF (LP.GT.0) CALL XERMSG ('SLATEC', 'LA05AS',
 
486
     *   'LENGTHS OF ARRAYS A(*) AND IND(*,2) ARE TOO SMALL.', -7, 1)
 
487
      G = -7.
 
488
      RETURN
 
489
C
 
490
  720 IPV = IPV + 1
 
491
      IW(IPV,1) = IR
 
492
      DO 730 I=1,N
 
493
         II = -IW(I,L+4)
 
494
         IF (II.GT.0) IW(II,1) = I
 
495
  730 CONTINUE
 
496
C
 
497
      IF (LP.GT.0) THEN
 
498
         XERN1 = 'ROWS'
 
499
         IF (L.EQ.2) XERN1 = 'COLUMNS'
 
500
         CALL XERMSG ('SLATEC', 'LA05AS', 'DEPENDANT ' // XERN1, -5, 1)
 
501
C
 
502
  740    WRITE (XERN1, '(I8)') IW(I,1)
 
503
         XERN2 = ' '
 
504
         IF (I+1.LE.IPV) WRITE (XERN2, '(I8)') IW(I+1,1)
 
505
         CALL XERMSG ('SLATEC', 'LA05AS',
 
506
     *      'DEPENDENT VECTOR INDICES ARE ' // XERN1 // ' AND ' //
 
507
     *      XERN2, -5, 1)
 
508
         I = I + 2
 
509
         IF (I.LE.IPV) GO TO 740
 
510
      ENDIF
 
511
      G = -5.
 
512
      RETURN
 
513
      END