001:       SUBROUTINE SLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
002:      $                   CNORM, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          DIAG, NORMIN, TRANS, UPLO
010:       INTEGER            INFO, N
011:       REAL               SCALE
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               AP( * ), CNORM( * ), X( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  SLATPS solves one of the triangular systems
021: *
022: *     A *x = s*b  or  A'*x = s*b
023: *
024: *  with scaling to prevent overflow, where A is an upper or lower
025: *  triangular matrix stored in packed form.  Here A' denotes the
026: *  transpose of A, x and b are n-element vectors, and s is a scaling
027: *  factor, usually less than or equal to 1, chosen so that the
028: *  components of x will be less than the overflow threshold.  If the
029: *  unscaled problem will not cause overflow, the Level 2 BLAS routine
030: *  STPSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
031: *  then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
032: *
033: *  Arguments
034: *  =========
035: *
036: *  UPLO    (input) CHARACTER*1
037: *          Specifies whether the matrix A is upper or lower triangular.
038: *          = 'U':  Upper triangular
039: *          = 'L':  Lower triangular
040: *
041: *  TRANS   (input) CHARACTER*1
042: *          Specifies the operation applied to A.
043: *          = 'N':  Solve A * x = s*b  (No transpose)
044: *          = 'T':  Solve A'* x = s*b  (Transpose)
045: *          = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)
046: *
047: *  DIAG    (input) CHARACTER*1
048: *          Specifies whether or not the matrix A is unit triangular.
049: *          = 'N':  Non-unit triangular
050: *          = 'U':  Unit triangular
051: *
052: *  NORMIN  (input) CHARACTER*1
053: *          Specifies whether CNORM has been set or not.
054: *          = 'Y':  CNORM contains the column norms on entry
055: *          = 'N':  CNORM is not set on entry.  On exit, the norms will
056: *                  be computed and stored in CNORM.
057: *
058: *  N       (input) INTEGER
059: *          The order of the matrix A.  N >= 0.
060: *
061: *  AP      (input) REAL array, dimension (N*(N+1)/2)
062: *          The upper or lower triangular matrix A, packed columnwise in
063: *          a linear array.  The j-th column of A is stored in the array
064: *          AP as follows:
065: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
066: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
067: *
068: *  X       (input/output) REAL array, dimension (N)
069: *          On entry, the right hand side b of the triangular system.
070: *          On exit, X is overwritten by the solution vector x.
071: *
072: *  SCALE   (output) REAL
073: *          The scaling factor s for the triangular system
074: *             A * x = s*b  or  A'* x = s*b.
075: *          If SCALE = 0, the matrix A is singular or badly scaled, and
076: *          the vector x is an exact or approximate solution to A*x = 0.
077: *
078: *  CNORM   (input or output) REAL array, dimension (N)
079: *
080: *          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
081: *          contains the norm of the off-diagonal part of the j-th column
082: *          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
083: *          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
084: *          must be greater than or equal to the 1-norm.
085: *
086: *          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
087: *          returns the 1-norm of the offdiagonal part of the j-th column
088: *          of A.
089: *
090: *  INFO    (output) INTEGER
091: *          = 0:  successful exit
092: *          < 0:  if INFO = -k, the k-th argument had an illegal value
093: *
094: *  Further Details
095: *  ======= =======
096: *
097: *  A rough bound on x is computed; if that is less than overflow, STPSV
098: *  is called, otherwise, specific code is used which checks for possible
099: *  overflow or divide-by-zero at every operation.
100: *
101: *  A columnwise scheme is used for solving A*x = b.  The basic algorithm
102: *  if A is lower triangular is
103: *
104: *       x[1:n] := b[1:n]
105: *       for j = 1, ..., n
106: *            x(j) := x(j) / A(j,j)
107: *            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
108: *       end
109: *
110: *  Define bounds on the components of x after j iterations of the loop:
111: *     M(j) = bound on x[1:j]
112: *     G(j) = bound on x[j+1:n]
113: *  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
114: *
115: *  Then for iteration j+1 we have
116: *     M(j+1) <= G(j) / | A(j+1,j+1) |
117: *     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
118: *            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
119: *
120: *  where CNORM(j+1) is greater than or equal to the infinity-norm of
121: *  column j+1 of A, not counting the diagonal.  Hence
122: *
123: *     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
124: *                  1<=i<=j
125: *  and
126: *
127: *     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
128: *                                   1<=i< j
129: *
130: *  Since |x(j)| <= M(j), we use the Level 2 BLAS routine STPSV if the
131: *  reciprocal of the largest M(j), j=1,..,n, is larger than
132: *  max(underflow, 1/overflow).
133: *
134: *  The bound on x(j) is also used to determine when a step in the
135: *  columnwise method can be performed without fear of overflow.  If
136: *  the computed bound is greater than a large constant, x is scaled to
137: *  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
138: *  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
139: *
140: *  Similarly, a row-wise scheme is used to solve A'*x = b.  The basic
141: *  algorithm for A upper triangular is
142: *
143: *       for j = 1, ..., n
144: *            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
145: *       end
146: *
147: *  We simultaneously compute two bounds
148: *       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
149: *       M(j) = bound on x(i), 1<=i<=j
150: *
151: *  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
152: *  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
153: *  Then the bound on x(j) is
154: *
155: *       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
156: *
157: *            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
158: *                      1<=i<=j
159: *
160: *  and we can safely call STPSV if 1/M(n) and 1/G(n) are both greater
161: *  than max(underflow, 1/overflow).
162: *
163: *  =====================================================================
164: *
165: *     .. Parameters ..
166:       REAL               ZERO, HALF, ONE
167:       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
168: *     ..
169: *     .. Local Scalars ..
170:       LOGICAL            NOTRAN, NOUNIT, UPPER
171:       INTEGER            I, IMAX, IP, J, JFIRST, JINC, JLAST, JLEN
172:       REAL               BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
173:      $                   TMAX, TSCAL, USCAL, XBND, XJ, XMAX
174: *     ..
175: *     .. External Functions ..
176:       LOGICAL            LSAME
177:       INTEGER            ISAMAX
178:       REAL               SASUM, SDOT, SLAMCH
179:       EXTERNAL           LSAME, ISAMAX, SASUM, SDOT, SLAMCH
180: *     ..
181: *     .. External Subroutines ..
182:       EXTERNAL           SAXPY, SSCAL, STPSV, XERBLA
183: *     ..
184: *     .. Intrinsic Functions ..
185:       INTRINSIC          ABS, MAX, MIN
186: *     ..
187: *     .. Executable Statements ..
188: *
189:       INFO = 0
190:       UPPER = LSAME( UPLO, 'U' )
191:       NOTRAN = LSAME( TRANS, 'N' )
192:       NOUNIT = LSAME( DIAG, 'N' )
193: *
194: *     Test the input parameters.
195: *
196:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
197:          INFO = -1
198:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
199:      $         LSAME( TRANS, 'C' ) ) THEN
200:          INFO = -2
201:       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
202:          INFO = -3
203:       ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
204:      $         LSAME( NORMIN, 'N' ) ) THEN
205:          INFO = -4
206:       ELSE IF( N.LT.0 ) THEN
207:          INFO = -5
208:       END IF
209:       IF( INFO.NE.0 ) THEN
210:          CALL XERBLA( 'SLATPS', -INFO )
211:          RETURN
212:       END IF
213: *
214: *     Quick return if possible
215: *
216:       IF( N.EQ.0 )
217:      $   RETURN
218: *
219: *     Determine machine dependent parameters to control overflow.
220: *
221:       SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )
222:       BIGNUM = ONE / SMLNUM
223:       SCALE = ONE
224: *
225:       IF( LSAME( NORMIN, 'N' ) ) THEN
226: *
227: *        Compute the 1-norm of each column, not including the diagonal.
228: *
229:          IF( UPPER ) THEN
230: *
231: *           A is upper triangular.
232: *
233:             IP = 1
234:             DO 10 J = 1, N
235:                CNORM( J ) = SASUM( J-1, AP( IP ), 1 )
236:                IP = IP + J
237:    10       CONTINUE
238:          ELSE
239: *
240: *           A is lower triangular.
241: *
242:             IP = 1
243:             DO 20 J = 1, N - 1
244:                CNORM( J ) = SASUM( N-J, AP( IP+1 ), 1 )
245:                IP = IP + N - J + 1
246:    20       CONTINUE
247:             CNORM( N ) = ZERO
248:          END IF
249:       END IF
250: *
251: *     Scale the column norms by TSCAL if the maximum element in CNORM is
252: *     greater than BIGNUM.
253: *
254:       IMAX = ISAMAX( N, CNORM, 1 )
255:       TMAX = CNORM( IMAX )
256:       IF( TMAX.LE.BIGNUM ) THEN
257:          TSCAL = ONE
258:       ELSE
259:          TSCAL = ONE / ( SMLNUM*TMAX )
260:          CALL SSCAL( N, TSCAL, CNORM, 1 )
261:       END IF
262: *
263: *     Compute a bound on the computed solution vector to see if the
264: *     Level 2 BLAS routine STPSV can be used.
265: *
266:       J = ISAMAX( N, X, 1 )
267:       XMAX = ABS( X( J ) )
268:       XBND = XMAX
269:       IF( NOTRAN ) THEN
270: *
271: *        Compute the growth in A * x = b.
272: *
273:          IF( UPPER ) THEN
274:             JFIRST = N
275:             JLAST = 1
276:             JINC = -1
277:          ELSE
278:             JFIRST = 1
279:             JLAST = N
280:             JINC = 1
281:          END IF
282: *
283:          IF( TSCAL.NE.ONE ) THEN
284:             GROW = ZERO
285:             GO TO 50
286:          END IF
287: *
288:          IF( NOUNIT ) THEN
289: *
290: *           A is non-unit triangular.
291: *
292: *           Compute GROW = 1/G(j) and XBND = 1/M(j).
293: *           Initially, G(0) = max{x(i), i=1,...,n}.
294: *
295:             GROW = ONE / MAX( XBND, SMLNUM )
296:             XBND = GROW
297:             IP = JFIRST*( JFIRST+1 ) / 2
298:             JLEN = N
299:             DO 30 J = JFIRST, JLAST, JINC
300: *
301: *              Exit the loop if the growth factor is too small.
302: *
303:                IF( GROW.LE.SMLNUM )
304:      $            GO TO 50
305: *
306: *              M(j) = G(j-1) / abs(A(j,j))
307: *
308:                TJJ = ABS( AP( IP ) )
309:                XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
310:                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
311: *
312: *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
313: *
314:                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
315:                ELSE
316: *
317: *                 G(j) could overflow, set GROW to 0.
318: *
319:                   GROW = ZERO
320:                END IF
321:                IP = IP + JINC*JLEN
322:                JLEN = JLEN - 1
323:    30       CONTINUE
324:             GROW = XBND
325:          ELSE
326: *
327: *           A is unit triangular.
328: *
329: *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
330: *
331:             GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
332:             DO 40 J = JFIRST, JLAST, JINC
333: *
334: *              Exit the loop if the growth factor is too small.
335: *
336:                IF( GROW.LE.SMLNUM )
337:      $            GO TO 50
338: *
339: *              G(j) = G(j-1)*( 1 + CNORM(j) )
340: *
341:                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
342:    40       CONTINUE
343:          END IF
344:    50    CONTINUE
345: *
346:       ELSE
347: *
348: *        Compute the growth in A' * x = b.
349: *
350:          IF( UPPER ) THEN
351:             JFIRST = 1
352:             JLAST = N
353:             JINC = 1
354:          ELSE
355:             JFIRST = N
356:             JLAST = 1
357:             JINC = -1
358:          END IF
359: *
360:          IF( TSCAL.NE.ONE ) THEN
361:             GROW = ZERO
362:             GO TO 80
363:          END IF
364: *
365:          IF( NOUNIT ) THEN
366: *
367: *           A is non-unit triangular.
368: *
369: *           Compute GROW = 1/G(j) and XBND = 1/M(j).
370: *           Initially, M(0) = max{x(i), i=1,...,n}.
371: *
372:             GROW = ONE / MAX( XBND, SMLNUM )
373:             XBND = GROW
374:             IP = JFIRST*( JFIRST+1 ) / 2
375:             JLEN = 1
376:             DO 60 J = JFIRST, JLAST, JINC
377: *
378: *              Exit the loop if the growth factor is too small.
379: *
380:                IF( GROW.LE.SMLNUM )
381:      $            GO TO 80
382: *
383: *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
384: *
385:                XJ = ONE + CNORM( J )
386:                GROW = MIN( GROW, XBND / XJ )
387: *
388: *              M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
389: *
390:                TJJ = ABS( AP( IP ) )
391:                IF( XJ.GT.TJJ )
392:      $            XBND = XBND*( TJJ / XJ )
393:                JLEN = JLEN + 1
394:                IP = IP + JINC*JLEN
395:    60       CONTINUE
396:             GROW = MIN( GROW, XBND )
397:          ELSE
398: *
399: *           A is unit triangular.
400: *
401: *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
402: *
403:             GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
404:             DO 70 J = JFIRST, JLAST, JINC
405: *
406: *              Exit the loop if the growth factor is too small.
407: *
408:                IF( GROW.LE.SMLNUM )
409:      $            GO TO 80
410: *
411: *              G(j) = ( 1 + CNORM(j) )*G(j-1)
412: *
413:                XJ = ONE + CNORM( J )
414:                GROW = GROW / XJ
415:    70       CONTINUE
416:          END IF
417:    80    CONTINUE
418:       END IF
419: *
420:       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
421: *
422: *        Use the Level 2 BLAS solve if the reciprocal of the bound on
423: *        elements of X is not too small.
424: *
425:          CALL STPSV( UPLO, TRANS, DIAG, N, AP, X, 1 )
426:       ELSE
427: *
428: *        Use a Level 1 BLAS solve, scaling intermediate results.
429: *
430:          IF( XMAX.GT.BIGNUM ) THEN
431: *
432: *           Scale X so that its components are less than or equal to
433: *           BIGNUM in absolute value.
434: *
435:             SCALE = BIGNUM / XMAX
436:             CALL SSCAL( N, SCALE, X, 1 )
437:             XMAX = BIGNUM
438:          END IF
439: *
440:          IF( NOTRAN ) THEN
441: *
442: *           Solve A * x = b
443: *
444:             IP = JFIRST*( JFIRST+1 ) / 2
445:             DO 100 J = JFIRST, JLAST, JINC
446: *
447: *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
448: *
449:                XJ = ABS( X( J ) )
450:                IF( NOUNIT ) THEN
451:                   TJJS = AP( IP )*TSCAL
452:                ELSE
453:                   TJJS = TSCAL
454:                   IF( TSCAL.EQ.ONE )
455:      $               GO TO 95
456:                END IF
457:                   TJJ = ABS( TJJS )
458:                   IF( TJJ.GT.SMLNUM ) THEN
459: *
460: *                    abs(A(j,j)) > SMLNUM:
461: *
462:                      IF( TJJ.LT.ONE ) THEN
463:                         IF( XJ.GT.TJJ*BIGNUM ) THEN
464: *
465: *                          Scale x by 1/b(j).
466: *
467:                            REC = ONE / XJ
468:                            CALL SSCAL( N, REC, X, 1 )
469:                            SCALE = SCALE*REC
470:                            XMAX = XMAX*REC
471:                         END IF
472:                      END IF
473:                      X( J ) = X( J ) / TJJS
474:                      XJ = ABS( X( J ) )
475:                   ELSE IF( TJJ.GT.ZERO ) THEN
476: *
477: *                    0 < abs(A(j,j)) <= SMLNUM:
478: *
479:                      IF( XJ.GT.TJJ*BIGNUM ) THEN
480: *
481: *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
482: *                       to avoid overflow when dividing by A(j,j).
483: *
484:                         REC = ( TJJ*BIGNUM ) / XJ
485:                         IF( CNORM( J ).GT.ONE ) THEN
486: *
487: *                          Scale by 1/CNORM(j) to avoid overflow when
488: *                          multiplying x(j) times column j.
489: *
490:                            REC = REC / CNORM( J )
491:                         END IF
492:                         CALL SSCAL( N, REC, X, 1 )
493:                         SCALE = SCALE*REC
494:                         XMAX = XMAX*REC
495:                      END IF
496:                      X( J ) = X( J ) / TJJS
497:                      XJ = ABS( X( J ) )
498:                   ELSE
499: *
500: *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
501: *                    scale = 0, and compute a solution to A*x = 0.
502: *
503:                      DO 90 I = 1, N
504:                         X( I ) = ZERO
505:    90                CONTINUE
506:                      X( J ) = ONE
507:                      XJ = ONE
508:                      SCALE = ZERO
509:                      XMAX = ZERO
510:                   END IF
511:    95          CONTINUE
512: *
513: *              Scale x if necessary to avoid overflow when adding a
514: *              multiple of column j of A.
515: *
516:                IF( XJ.GT.ONE ) THEN
517:                   REC = ONE / XJ
518:                   IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
519: *
520: *                    Scale x by 1/(2*abs(x(j))).
521: *
522:                      REC = REC*HALF
523:                      CALL SSCAL( N, REC, X, 1 )
524:                      SCALE = SCALE*REC
525:                   END IF
526:                ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
527: *
528: *                 Scale x by 1/2.
529: *
530:                   CALL SSCAL( N, HALF, X, 1 )
531:                   SCALE = SCALE*HALF
532:                END IF
533: *
534:                IF( UPPER ) THEN
535:                   IF( J.GT.1 ) THEN
536: *
537: *                    Compute the update
538: *                       x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
539: *
540:                      CALL SAXPY( J-1, -X( J )*TSCAL, AP( IP-J+1 ), 1, X,
541:      $                           1 )
542:                      I = ISAMAX( J-1, X, 1 )
543:                      XMAX = ABS( X( I ) )
544:                   END IF
545:                   IP = IP - J
546:                ELSE
547:                   IF( J.LT.N ) THEN
548: *
549: *                    Compute the update
550: *                       x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
551: *
552:                      CALL SAXPY( N-J, -X( J )*TSCAL, AP( IP+1 ), 1,
553:      $                           X( J+1 ), 1 )
554:                      I = J + ISAMAX( N-J, X( J+1 ), 1 )
555:                      XMAX = ABS( X( I ) )
556:                   END IF
557:                   IP = IP + N - J + 1
558:                END IF
559:   100       CONTINUE
560: *
561:          ELSE
562: *
563: *           Solve A' * x = b
564: *
565:             IP = JFIRST*( JFIRST+1 ) / 2
566:             JLEN = 1
567:             DO 140 J = JFIRST, JLAST, JINC
568: *
569: *              Compute x(j) = b(j) - sum A(k,j)*x(k).
570: *                                    k<>j
571: *
572:                XJ = ABS( X( J ) )
573:                USCAL = TSCAL
574:                REC = ONE / MAX( XMAX, ONE )
575:                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
576: *
577: *                 If x(j) could overflow, scale x by 1/(2*XMAX).
578: *
579:                   REC = REC*HALF
580:                   IF( NOUNIT ) THEN
581:                      TJJS = AP( IP )*TSCAL
582:                   ELSE
583:                      TJJS = TSCAL
584:                   END IF
585:                      TJJ = ABS( TJJS )
586:                      IF( TJJ.GT.ONE ) THEN
587: *
588: *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
589: *
590:                         REC = MIN( ONE, REC*TJJ )
591:                         USCAL = USCAL / TJJS
592:                      END IF
593:                   IF( REC.LT.ONE ) THEN
594:                      CALL SSCAL( N, REC, X, 1 )
595:                      SCALE = SCALE*REC
596:                      XMAX = XMAX*REC
597:                   END IF
598:                END IF
599: *
600:                SUMJ = ZERO
601:                IF( USCAL.EQ.ONE ) THEN
602: *
603: *                 If the scaling needed for A in the dot product is 1,
604: *                 call SDOT to perform the dot product.
605: *
606:                   IF( UPPER ) THEN
607:                      SUMJ = SDOT( J-1, AP( IP-J+1 ), 1, X, 1 )
608:                   ELSE IF( J.LT.N ) THEN
609:                      SUMJ = SDOT( N-J, AP( IP+1 ), 1, X( J+1 ), 1 )
610:                   END IF
611:                ELSE
612: *
613: *                 Otherwise, use in-line code for the dot product.
614: *
615:                   IF( UPPER ) THEN
616:                      DO 110 I = 1, J - 1
617:                         SUMJ = SUMJ + ( AP( IP-J+I )*USCAL )*X( I )
618:   110                CONTINUE
619:                   ELSE IF( J.LT.N ) THEN
620:                      DO 120 I = 1, N - J
621:                         SUMJ = SUMJ + ( AP( IP+I )*USCAL )*X( J+I )
622:   120                CONTINUE
623:                   END IF
624:                END IF
625: *
626:                IF( USCAL.EQ.TSCAL ) THEN
627: *
628: *                 Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
629: *                 was not used to scale the dotproduct.
630: *
631:                   X( J ) = X( J ) - SUMJ
632:                   XJ = ABS( X( J ) )
633:                   IF( NOUNIT ) THEN
634: *
635: *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
636: *
637:                      TJJS = AP( IP )*TSCAL
638:                   ELSE
639:                      TJJS = TSCAL
640:                      IF( TSCAL.EQ.ONE )
641:      $                  GO TO 135
642:                   END IF
643:                      TJJ = ABS( TJJS )
644:                      IF( TJJ.GT.SMLNUM ) THEN
645: *
646: *                       abs(A(j,j)) > SMLNUM:
647: *
648:                         IF( TJJ.LT.ONE ) THEN
649:                            IF( XJ.GT.TJJ*BIGNUM ) THEN
650: *
651: *                             Scale X by 1/abs(x(j)).
652: *
653:                               REC = ONE / XJ
654:                               CALL SSCAL( N, REC, X, 1 )
655:                               SCALE = SCALE*REC
656:                               XMAX = XMAX*REC
657:                            END IF
658:                         END IF
659:                         X( J ) = X( J ) / TJJS
660:                      ELSE IF( TJJ.GT.ZERO ) THEN
661: *
662: *                       0 < abs(A(j,j)) <= SMLNUM:
663: *
664:                         IF( XJ.GT.TJJ*BIGNUM ) THEN
665: *
666: *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
667: *
668:                            REC = ( TJJ*BIGNUM ) / XJ
669:                            CALL SSCAL( N, REC, X, 1 )
670:                            SCALE = SCALE*REC
671:                            XMAX = XMAX*REC
672:                         END IF
673:                         X( J ) = X( J ) / TJJS
674:                      ELSE
675: *
676: *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
677: *                       scale = 0, and compute a solution to A'*x = 0.
678: *
679:                         DO 130 I = 1, N
680:                            X( I ) = ZERO
681:   130                   CONTINUE
682:                         X( J ) = ONE
683:                         SCALE = ZERO
684:                         XMAX = ZERO
685:                      END IF
686:   135             CONTINUE
687:                ELSE
688: *
689: *                 Compute x(j) := x(j) / A(j,j)  - sumj if the dot
690: *                 product has already been divided by 1/A(j,j).
691: *
692:                   X( J ) = X( J ) / TJJS - SUMJ
693:                END IF
694:                XMAX = MAX( XMAX, ABS( X( J ) ) )
695:                JLEN = JLEN + 1
696:                IP = IP + JINC*JLEN
697:   140       CONTINUE
698:          END IF
699:          SCALE = SCALE / TSCAL
700:       END IF
701: *
702: *     Scale the column norms by 1/TSCAL for return.
703: *
704:       IF( TSCAL.NE.ONE ) THEN
705:          CALL SSCAL( N, ONE / TSCAL, CNORM, 1 )
706:       END IF
707: *
708:       RETURN
709: *
710: *     End of SLATPS
711: *
712:       END
713: