LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ clatps()

subroutine clatps ( character  uplo,
character  trans,
character  diag,
character  normin,
integer  n,
complex, dimension( * )  ap,
complex, dimension( * )  x,
real  scale,
real, dimension( * )  cnorm,
integer  info 
)

CLATPS solves a triangular system of equations with the matrix held in packed storage.

Download CLATPS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CLATPS solves one of the triangular systems

    A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,

 with scaling to prevent overflow, where A is an upper or lower
 triangular matrix stored in packed form.  Here A**T denotes the
 transpose of A, A**H denotes the conjugate transpose of A, x and b
 are n-element vectors, and s is a scaling factor, usually less than
 or equal to 1, chosen so that the components of x will be less than
 the overflow threshold.  If the unscaled problem will not cause
 overflow, the Level 2 BLAS routine CTPSV is called. If the matrix A
 is singular (A(j,j) = 0 for some j), then s is set to 0 and a
 non-trivial solution to A*x = 0 is returned.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the matrix A is upper or lower triangular.
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the operation applied to A.
          = 'N':  Solve A * x = s*b     (No transpose)
          = 'T':  Solve A**T * x = s*b  (Transpose)
          = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
[in]DIAG
          DIAG is CHARACTER*1
          Specifies whether or not the matrix A is unit triangular.
          = 'N':  Non-unit triangular
          = 'U':  Unit triangular
[in]NORMIN
          NORMIN is CHARACTER*1
          Specifies whether CNORM has been set or not.
          = 'Y':  CNORM contains the column norms on entry
          = 'N':  CNORM is not set on entry.  On exit, the norms will
                  be computed and stored in CNORM.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]AP
          AP is COMPLEX array, dimension (N*(N+1)/2)
          The upper or lower triangular matrix A, packed columnwise in
          a linear array.  The j-th column of A is stored in the array
          AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
[in,out]X
          X is COMPLEX array, dimension (N)
          On entry, the right hand side b of the triangular system.
          On exit, X is overwritten by the solution vector x.
[out]SCALE
          SCALE is REAL
          The scaling factor s for the triangular system
             A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b.
          If SCALE = 0, the matrix A is singular or badly scaled, and
          the vector x is an exact or approximate solution to A*x = 0.
[in,out]CNORM
          CNORM is REAL array, dimension (N)

          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
          contains the norm of the off-diagonal part of the j-th column
          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
          must be greater than or equal to the 1-norm.

          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
          returns the 1-norm of the offdiagonal part of the j-th column
          of A.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -k, the k-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  A rough bound on x is computed; if that is less than overflow, CTPSV
  is called, otherwise, specific code is used which checks for possible
  overflow or divide-by-zero at every operation.

  A columnwise scheme is used for solving A*x = b.  The basic algorithm
  if A is lower triangular is

       x[1:n] := b[1:n]
       for j = 1, ..., n
            x(j) := x(j) / A(j,j)
            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
       end

  Define bounds on the components of x after j iterations of the loop:
     M(j) = bound on x[1:j]
     G(j) = bound on x[j+1:n]
  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.

  Then for iteration j+1 we have
     M(j+1) <= G(j) / | A(j+1,j+1) |
     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )

  where CNORM(j+1) is greater than or equal to the infinity-norm of
  column j+1 of A, not counting the diagonal.  Hence

     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
                  1<=i<=j
  and

     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
                                   1<=i< j

  Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTPSV if the
  reciprocal of the largest M(j), j=1,..,n, is larger than
  max(underflow, 1/overflow).

  The bound on x(j) is also used to determine when a step in the
  columnwise method can be performed without fear of overflow.  If
  the computed bound is greater than a large constant, x is scaled to
  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.

  Similarly, a row-wise scheme is used to solve A**T *x = b  or
  A**H *x = b.  The basic algorithm for A upper triangular is

       for j = 1, ..., n
            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
       end

  We simultaneously compute two bounds
       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
       M(j) = bound on x(i), 1<=i<=j

  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
  Then the bound on x(j) is

       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |

            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
                      1<=i<=j

  and we can safely call CTPSV if 1/M(n) and 1/G(n) are both greater
  than max(underflow, 1/overflow).

Definition at line 229 of file clatps.f.

231*
232* -- LAPACK auxiliary routine --
233* -- LAPACK is a software package provided by Univ. of Tennessee, --
234* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
235*
236* .. Scalar Arguments ..
237 CHARACTER DIAG, NORMIN, TRANS, UPLO
238 INTEGER INFO, N
239 REAL SCALE
240* ..
241* .. Array Arguments ..
242 REAL CNORM( * )
243 COMPLEX AP( * ), X( * )
244* ..
245*
246* =====================================================================
247*
248* .. Parameters ..
249 REAL ZERO, HALF, ONE, TWO
250 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
251 $ two = 2.0e+0 )
252* ..
253* .. Local Scalars ..
254 LOGICAL NOTRAN, NOUNIT, UPPER
255 INTEGER I, IMAX, IP, J, JFIRST, JINC, JLAST, JLEN
256 REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
257 $ XBND, XJ, XMAX
258 COMPLEX CSUMJ, TJJS, USCAL, ZDUM
259* ..
260* .. External Functions ..
261 LOGICAL LSAME
262 INTEGER ICAMAX, ISAMAX
263 REAL SCASUM, SLAMCH
264 COMPLEX CDOTC, CDOTU, CLADIV
265 EXTERNAL lsame, icamax, isamax, scasum, slamch, cdotc,
266 $ cdotu, cladiv
267* ..
268* .. External Subroutines ..
269 EXTERNAL caxpy, csscal, ctpsv, sscal, xerbla
270* ..
271* .. Intrinsic Functions ..
272 INTRINSIC abs, aimag, cmplx, conjg, max, min, real
273* ..
274* .. Statement Functions ..
275 REAL CABS1, CABS2
276* ..
277* .. Statement Function definitions ..
278 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
279 cabs2( zdum ) = abs( real( zdum ) / 2. ) +
280 $ abs( aimag( zdum ) / 2. )
281* ..
282* .. Executable Statements ..
283*
284 info = 0
285 upper = lsame( uplo, 'U' )
286 notran = lsame( trans, 'N' )
287 nounit = lsame( diag, 'N' )
288*
289* Test the input parameters.
290*
291 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
292 info = -1
293 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
294 $ lsame( trans, 'C' ) ) THEN
295 info = -2
296 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
297 info = -3
298 ELSE IF( .NOT.lsame( normin, 'Y' ) .AND. .NOT.
299 $ lsame( normin, 'N' ) ) THEN
300 info = -4
301 ELSE IF( n.LT.0 ) THEN
302 info = -5
303 END IF
304 IF( info.NE.0 ) THEN
305 CALL xerbla( 'CLATPS', -info )
306 RETURN
307 END IF
308*
309* Quick return if possible
310*
311 IF( n.EQ.0 )
312 $ RETURN
313*
314* Determine machine dependent parameters to control overflow.
315*
316 smlnum = slamch( 'Safe minimum' )
317 bignum = one / smlnum
318 smlnum = smlnum / slamch( 'Precision' )
319 bignum = one / smlnum
320 scale = one
321*
322 IF( lsame( normin, 'N' ) ) THEN
323*
324* Compute the 1-norm of each column, not including the diagonal.
325*
326 IF( upper ) THEN
327*
328* A is upper triangular.
329*
330 ip = 1
331 DO 10 j = 1, n
332 cnorm( j ) = scasum( j-1, ap( ip ), 1 )
333 ip = ip + j
334 10 CONTINUE
335 ELSE
336*
337* A is lower triangular.
338*
339 ip = 1
340 DO 20 j = 1, n - 1
341 cnorm( j ) = scasum( n-j, ap( ip+1 ), 1 )
342 ip = ip + n - j + 1
343 20 CONTINUE
344 cnorm( n ) = zero
345 END IF
346 END IF
347*
348* Scale the column norms by TSCAL if the maximum element in CNORM is
349* greater than BIGNUM/2.
350*
351 imax = isamax( n, cnorm, 1 )
352 tmax = cnorm( imax )
353 IF( tmax.LE.bignum*half ) THEN
354 tscal = one
355 ELSE
356 tscal = half / ( smlnum*tmax )
357 CALL sscal( n, tscal, cnorm, 1 )
358 END IF
359*
360* Compute a bound on the computed solution vector to see if the
361* Level 2 BLAS routine CTPSV can be used.
362*
363 xmax = zero
364 DO 30 j = 1, n
365 xmax = max( xmax, cabs2( x( j ) ) )
366 30 CONTINUE
367 xbnd = xmax
368 IF( notran ) THEN
369*
370* Compute the growth in A * x = b.
371*
372 IF( upper ) THEN
373 jfirst = n
374 jlast = 1
375 jinc = -1
376 ELSE
377 jfirst = 1
378 jlast = n
379 jinc = 1
380 END IF
381*
382 IF( tscal.NE.one ) THEN
383 grow = zero
384 GO TO 60
385 END IF
386*
387 IF( nounit ) THEN
388*
389* A is non-unit triangular.
390*
391* Compute GROW = 1/G(j) and XBND = 1/M(j).
392* Initially, G(0) = max{x(i), i=1,...,n}.
393*
394 grow = half / max( xbnd, smlnum )
395 xbnd = grow
396 ip = jfirst*( jfirst+1 ) / 2
397 jlen = n
398 DO 40 j = jfirst, jlast, jinc
399*
400* Exit the loop if the growth factor is too small.
401*
402 IF( grow.LE.smlnum )
403 $ GO TO 60
404*
405 tjjs = ap( ip )
406 tjj = cabs1( tjjs )
407*
408 IF( tjj.GE.smlnum ) THEN
409*
410* M(j) = G(j-1) / abs(A(j,j))
411*
412 xbnd = min( xbnd, min( one, tjj )*grow )
413 ELSE
414*
415* M(j) could overflow, set XBND to 0.
416*
417 xbnd = zero
418 END IF
419*
420 IF( tjj+cnorm( j ).GE.smlnum ) THEN
421*
422* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
423*
424 grow = grow*( tjj / ( tjj+cnorm( j ) ) )
425 ELSE
426*
427* G(j) could overflow, set GROW to 0.
428*
429 grow = zero
430 END IF
431 ip = ip + jinc*jlen
432 jlen = jlen - 1
433 40 CONTINUE
434 grow = xbnd
435 ELSE
436*
437* A is unit triangular.
438*
439* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
440*
441 grow = min( one, half / max( xbnd, smlnum ) )
442 DO 50 j = jfirst, jlast, jinc
443*
444* Exit the loop if the growth factor is too small.
445*
446 IF( grow.LE.smlnum )
447 $ GO TO 60
448*
449* G(j) = G(j-1)*( 1 + CNORM(j) )
450*
451 grow = grow*( one / ( one+cnorm( j ) ) )
452 50 CONTINUE
453 END IF
454 60 CONTINUE
455*
456 ELSE
457*
458* Compute the growth in A**T * x = b or A**H * x = b.
459*
460 IF( upper ) THEN
461 jfirst = 1
462 jlast = n
463 jinc = 1
464 ELSE
465 jfirst = n
466 jlast = 1
467 jinc = -1
468 END IF
469*
470 IF( tscal.NE.one ) THEN
471 grow = zero
472 GO TO 90
473 END IF
474*
475 IF( nounit ) THEN
476*
477* A is non-unit triangular.
478*
479* Compute GROW = 1/G(j) and XBND = 1/M(j).
480* Initially, M(0) = max{x(i), i=1,...,n}.
481*
482 grow = half / max( xbnd, smlnum )
483 xbnd = grow
484 ip = jfirst*( jfirst+1 ) / 2
485 jlen = 1
486 DO 70 j = jfirst, jlast, jinc
487*
488* Exit the loop if the growth factor is too small.
489*
490 IF( grow.LE.smlnum )
491 $ GO TO 90
492*
493* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
494*
495 xj = one + cnorm( j )
496 grow = min( grow, xbnd / xj )
497*
498 tjjs = ap( ip )
499 tjj = cabs1( tjjs )
500*
501 IF( tjj.GE.smlnum ) THEN
502*
503* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
504*
505 IF( xj.GT.tjj )
506 $ xbnd = xbnd*( tjj / xj )
507 ELSE
508*
509* M(j) could overflow, set XBND to 0.
510*
511 xbnd = zero
512 END IF
513 jlen = jlen + 1
514 ip = ip + jinc*jlen
515 70 CONTINUE
516 grow = min( grow, xbnd )
517 ELSE
518*
519* A is unit triangular.
520*
521* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
522*
523 grow = min( one, half / max( xbnd, smlnum ) )
524 DO 80 j = jfirst, jlast, jinc
525*
526* Exit the loop if the growth factor is too small.
527*
528 IF( grow.LE.smlnum )
529 $ GO TO 90
530*
531* G(j) = ( 1 + CNORM(j) )*G(j-1)
532*
533 xj = one + cnorm( j )
534 grow = grow / xj
535 80 CONTINUE
536 END IF
537 90 CONTINUE
538 END IF
539*
540 IF( ( grow*tscal ).GT.smlnum ) THEN
541*
542* Use the Level 2 BLAS solve if the reciprocal of the bound on
543* elements of X is not too small.
544*
545 CALL ctpsv( uplo, trans, diag, n, ap, x, 1 )
546 ELSE
547*
548* Use a Level 1 BLAS solve, scaling intermediate results.
549*
550 IF( xmax.GT.bignum*half ) THEN
551*
552* Scale X so that its components are less than or equal to
553* BIGNUM in absolute value.
554*
555 scale = ( bignum*half ) / xmax
556 CALL csscal( n, scale, x, 1 )
557 xmax = bignum
558 ELSE
559 xmax = xmax*two
560 END IF
561*
562 IF( notran ) THEN
563*
564* Solve A * x = b
565*
566 ip = jfirst*( jfirst+1 ) / 2
567 DO 110 j = jfirst, jlast, jinc
568*
569* Compute x(j) = b(j) / A(j,j), scaling x if necessary.
570*
571 xj = cabs1( x( j ) )
572 IF( nounit ) THEN
573 tjjs = ap( ip )*tscal
574 ELSE
575 tjjs = tscal
576 IF( tscal.EQ.one )
577 $ GO TO 105
578 END IF
579 tjj = cabs1( tjjs )
580 IF( tjj.GT.smlnum ) THEN
581*
582* abs(A(j,j)) > SMLNUM:
583*
584 IF( tjj.LT.one ) THEN
585 IF( xj.GT.tjj*bignum ) THEN
586*
587* Scale x by 1/b(j).
588*
589 rec = one / xj
590 CALL csscal( n, rec, x, 1 )
591 scale = scale*rec
592 xmax = xmax*rec
593 END IF
594 END IF
595 x( j ) = cladiv( x( j ), tjjs )
596 xj = cabs1( x( j ) )
597 ELSE IF( tjj.GT.zero ) THEN
598*
599* 0 < abs(A(j,j)) <= SMLNUM:
600*
601 IF( xj.GT.tjj*bignum ) THEN
602*
603* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
604* to avoid overflow when dividing by A(j,j).
605*
606 rec = ( tjj*bignum ) / xj
607 IF( cnorm( j ).GT.one ) THEN
608*
609* Scale by 1/CNORM(j) to avoid overflow when
610* multiplying x(j) times column j.
611*
612 rec = rec / cnorm( j )
613 END IF
614 CALL csscal( n, rec, x, 1 )
615 scale = scale*rec
616 xmax = xmax*rec
617 END IF
618 x( j ) = cladiv( x( j ), tjjs )
619 xj = cabs1( x( j ) )
620 ELSE
621*
622* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
623* scale = 0, and compute a solution to A*x = 0.
624*
625 DO 100 i = 1, n
626 x( i ) = zero
627 100 CONTINUE
628 x( j ) = one
629 xj = one
630 scale = zero
631 xmax = zero
632 END IF
633 105 CONTINUE
634*
635* Scale x if necessary to avoid overflow when adding a
636* multiple of column j of A.
637*
638 IF( xj.GT.one ) THEN
639 rec = one / xj
640 IF( cnorm( j ).GT.( bignum-xmax )*rec ) THEN
641*
642* Scale x by 1/(2*abs(x(j))).
643*
644 rec = rec*half
645 CALL csscal( n, rec, x, 1 )
646 scale = scale*rec
647 END IF
648 ELSE IF( xj*cnorm( j ).GT.( bignum-xmax ) ) THEN
649*
650* Scale x by 1/2.
651*
652 CALL csscal( n, half, x, 1 )
653 scale = scale*half
654 END IF
655*
656 IF( upper ) THEN
657 IF( j.GT.1 ) THEN
658*
659* Compute the update
660* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
661*
662 CALL caxpy( j-1, -x( j )*tscal, ap( ip-j+1 ), 1, x,
663 $ 1 )
664 i = icamax( j-1, x, 1 )
665 xmax = cabs1( x( i ) )
666 END IF
667 ip = ip - j
668 ELSE
669 IF( j.LT.n ) THEN
670*
671* Compute the update
672* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
673*
674 CALL caxpy( n-j, -x( j )*tscal, ap( ip+1 ), 1,
675 $ x( j+1 ), 1 )
676 i = j + icamax( n-j, x( j+1 ), 1 )
677 xmax = cabs1( x( i ) )
678 END IF
679 ip = ip + n - j + 1
680 END IF
681 110 CONTINUE
682*
683 ELSE IF( lsame( trans, 'T' ) ) THEN
684*
685* Solve A**T * x = b
686*
687 ip = jfirst*( jfirst+1 ) / 2
688 jlen = 1
689 DO 150 j = jfirst, jlast, jinc
690*
691* Compute x(j) = b(j) - sum A(k,j)*x(k).
692* k<>j
693*
694 xj = cabs1( x( j ) )
695 uscal = tscal
696 rec = one / max( xmax, one )
697 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
698*
699* If x(j) could overflow, scale x by 1/(2*XMAX).
700*
701 rec = rec*half
702 IF( nounit ) THEN
703 tjjs = ap( ip )*tscal
704 ELSE
705 tjjs = tscal
706 END IF
707 tjj = cabs1( tjjs )
708 IF( tjj.GT.one ) THEN
709*
710* Divide by A(j,j) when scaling x if A(j,j) > 1.
711*
712 rec = min( one, rec*tjj )
713 uscal = cladiv( uscal, tjjs )
714 END IF
715 IF( rec.LT.one ) THEN
716 CALL csscal( n, rec, x, 1 )
717 scale = scale*rec
718 xmax = xmax*rec
719 END IF
720 END IF
721*
722 csumj = zero
723 IF( uscal.EQ.cmplx( one ) ) THEN
724*
725* If the scaling needed for A in the dot product is 1,
726* call CDOTU to perform the dot product.
727*
728 IF( upper ) THEN
729 csumj = cdotu( j-1, ap( ip-j+1 ), 1, x, 1 )
730 ELSE IF( j.LT.n ) THEN
731 csumj = cdotu( n-j, ap( ip+1 ), 1, x( j+1 ), 1 )
732 END IF
733 ELSE
734*
735* Otherwise, use in-line code for the dot product.
736*
737 IF( upper ) THEN
738 DO 120 i = 1, j - 1
739 csumj = csumj + ( ap( ip-j+i )*uscal )*x( i )
740 120 CONTINUE
741 ELSE IF( j.LT.n ) THEN
742 DO 130 i = 1, n - j
743 csumj = csumj + ( ap( ip+i )*uscal )*x( j+i )
744 130 CONTINUE
745 END IF
746 END IF
747*
748 IF( uscal.EQ.cmplx( tscal ) ) THEN
749*
750* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
751* was not used to scale the dotproduct.
752*
753 x( j ) = x( j ) - csumj
754 xj = cabs1( x( j ) )
755 IF( nounit ) THEN
756*
757* Compute x(j) = x(j) / A(j,j), scaling if necessary.
758*
759 tjjs = ap( ip )*tscal
760 ELSE
761 tjjs = tscal
762 IF( tscal.EQ.one )
763 $ GO TO 145
764 END IF
765 tjj = cabs1( tjjs )
766 IF( tjj.GT.smlnum ) THEN
767*
768* abs(A(j,j)) > SMLNUM:
769*
770 IF( tjj.LT.one ) THEN
771 IF( xj.GT.tjj*bignum ) THEN
772*
773* Scale X by 1/abs(x(j)).
774*
775 rec = one / xj
776 CALL csscal( n, rec, x, 1 )
777 scale = scale*rec
778 xmax = xmax*rec
779 END IF
780 END IF
781 x( j ) = cladiv( x( j ), tjjs )
782 ELSE IF( tjj.GT.zero ) THEN
783*
784* 0 < abs(A(j,j)) <= SMLNUM:
785*
786 IF( xj.GT.tjj*bignum ) THEN
787*
788* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
789*
790 rec = ( tjj*bignum ) / xj
791 CALL csscal( n, rec, x, 1 )
792 scale = scale*rec
793 xmax = xmax*rec
794 END IF
795 x( j ) = cladiv( x( j ), tjjs )
796 ELSE
797*
798* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
799* scale = 0 and compute a solution to A**T *x = 0.
800*
801 DO 140 i = 1, n
802 x( i ) = zero
803 140 CONTINUE
804 x( j ) = one
805 scale = zero
806 xmax = zero
807 END IF
808 145 CONTINUE
809 ELSE
810*
811* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
812* product has already been divided by 1/A(j,j).
813*
814 x( j ) = cladiv( x( j ), tjjs ) - csumj
815 END IF
816 xmax = max( xmax, cabs1( x( j ) ) )
817 jlen = jlen + 1
818 ip = ip + jinc*jlen
819 150 CONTINUE
820*
821 ELSE
822*
823* Solve A**H * x = b
824*
825 ip = jfirst*( jfirst+1 ) / 2
826 jlen = 1
827 DO 190 j = jfirst, jlast, jinc
828*
829* Compute x(j) = b(j) - sum A(k,j)*x(k).
830* k<>j
831*
832 xj = cabs1( x( j ) )
833 uscal = tscal
834 rec = one / max( xmax, one )
835 IF( cnorm( j ).GT.( bignum-xj )*rec ) THEN
836*
837* If x(j) could overflow, scale x by 1/(2*XMAX).
838*
839 rec = rec*half
840 IF( nounit ) THEN
841 tjjs = conjg( ap( ip ) )*tscal
842 ELSE
843 tjjs = tscal
844 END IF
845 tjj = cabs1( tjjs )
846 IF( tjj.GT.one ) THEN
847*
848* Divide by A(j,j) when scaling x if A(j,j) > 1.
849*
850 rec = min( one, rec*tjj )
851 uscal = cladiv( uscal, tjjs )
852 END IF
853 IF( rec.LT.one ) THEN
854 CALL csscal( n, rec, x, 1 )
855 scale = scale*rec
856 xmax = xmax*rec
857 END IF
858 END IF
859*
860 csumj = zero
861 IF( uscal.EQ.cmplx( one ) ) THEN
862*
863* If the scaling needed for A in the dot product is 1,
864* call CDOTC to perform the dot product.
865*
866 IF( upper ) THEN
867 csumj = cdotc( j-1, ap( ip-j+1 ), 1, x, 1 )
868 ELSE IF( j.LT.n ) THEN
869 csumj = cdotc( n-j, ap( ip+1 ), 1, x( j+1 ), 1 )
870 END IF
871 ELSE
872*
873* Otherwise, use in-line code for the dot product.
874*
875 IF( upper ) THEN
876 DO 160 i = 1, j - 1
877 csumj = csumj + ( conjg( ap( ip-j+i ) )*uscal )*
878 $ x( i )
879 160 CONTINUE
880 ELSE IF( j.LT.n ) THEN
881 DO 170 i = 1, n - j
882 csumj = csumj + ( conjg( ap( ip+i ) )*uscal )*
883 $ x( j+i )
884 170 CONTINUE
885 END IF
886 END IF
887*
888 IF( uscal.EQ.cmplx( tscal ) ) THEN
889*
890* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
891* was not used to scale the dotproduct.
892*
893 x( j ) = x( j ) - csumj
894 xj = cabs1( x( j ) )
895 IF( nounit ) THEN
896*
897* Compute x(j) = x(j) / A(j,j), scaling if necessary.
898*
899 tjjs = conjg( ap( ip ) )*tscal
900 ELSE
901 tjjs = tscal
902 IF( tscal.EQ.one )
903 $ GO TO 185
904 END IF
905 tjj = cabs1( tjjs )
906 IF( tjj.GT.smlnum ) THEN
907*
908* abs(A(j,j)) > SMLNUM:
909*
910 IF( tjj.LT.one ) THEN
911 IF( xj.GT.tjj*bignum ) THEN
912*
913* Scale X by 1/abs(x(j)).
914*
915 rec = one / xj
916 CALL csscal( n, rec, x, 1 )
917 scale = scale*rec
918 xmax = xmax*rec
919 END IF
920 END IF
921 x( j ) = cladiv( x( j ), tjjs )
922 ELSE IF( tjj.GT.zero ) THEN
923*
924* 0 < abs(A(j,j)) <= SMLNUM:
925*
926 IF( xj.GT.tjj*bignum ) THEN
927*
928* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
929*
930 rec = ( tjj*bignum ) / xj
931 CALL csscal( n, rec, x, 1 )
932 scale = scale*rec
933 xmax = xmax*rec
934 END IF
935 x( j ) = cladiv( x( j ), tjjs )
936 ELSE
937*
938* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
939* scale = 0 and compute a solution to A**H *x = 0.
940*
941 DO 180 i = 1, n
942 x( i ) = zero
943 180 CONTINUE
944 x( j ) = one
945 scale = zero
946 xmax = zero
947 END IF
948 185 CONTINUE
949 ELSE
950*
951* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
952* product has already been divided by 1/A(j,j).
953*
954 x( j ) = cladiv( x( j ), tjjs ) - csumj
955 END IF
956 xmax = max( xmax, cabs1( x( j ) ) )
957 jlen = jlen + 1
958 ip = ip + jinc*jlen
959 190 CONTINUE
960 END IF
961 scale = scale / tscal
962 END IF
963*
964* Scale the column norms by 1/TSCAL for return.
965*
966 IF( tscal.NE.one ) THEN
967 CALL sscal( n, one / tscal, cnorm, 1 )
968 END IF
969*
970 RETURN
971*
972* End of CLATPS
973*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
real function scasum(n, cx, incx)
SCASUM
Definition scasum.f:72
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
complex function cdotu(n, cx, incx, cy, incy)
CDOTU
Definition cdotu.f:83
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
complex function cladiv(x, y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition cladiv.f:64
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine ctpsv(uplo, trans, diag, n, ap, x, incx)
CTPSV
Definition ctpsv.f:144
Here is the call graph for this function:
Here is the caller graph for this function: