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