src/SOLIDE/solid.f

   1       subroutine transit(compo,dt0)
   2       include "calcium.hf"
   3       double precision a(5,24),b(24),tn(24)
   4       double precision q(6)
   5       dimension maille(2),t(24),puiss(24)
   6       dimension text(6),rflu(6)
   7       dimension tpar(6),rpar(6)
   8       double precision dt0
   9       integer compo
  10
  11       dt=dt0
  12
  13       ldab=5
  14       ldb=24
  15       nrhs=1
  16       npt=0
  17       r=1.
  18       ro=1.
  19       rext=0.5
  20       te=10.
  21       do i=1,24
  22       tn(i)=100.
  23       t(i)=100.
  24       puiss(i)=100.
  25       enddo
  26       do j=1,6
  27       q(j)=50.
  28       enddo
  29 c
  30 c   construction de la matrice (laplacien)
  31 c
  32       do i = 1, 5
  33         do j = 1, 24
  34           a(i,j) = 0.
  35         enddo
  36       enddo
  37
  38       do i=2,3
  39       do j=2,5
  40       npt=i+(j-1)*4
  41       a(1,npt)=4./r+ro/dt
  42       a(2,npt)=-1./r
  43       a(5,npt)=-1./r
  44       enddo
  45       enddo
  46       do i=2,3
  47       npt=i
  48       a(1,npt)=3./r+ro/dt
  49       a(2,npt)=-1./r
  50       a(5,npt)=-1./r
  51       npt=i+20
  52       a(1,npt)=3./r+ro/dt
  53       a(2,npt)=-1./r
  54       a(5,npt)=0.
  55       enddo
  56       do j=2,5
  57       npt=1+(j-1)*4
  58       a(1,npt)=3./r+ro/dt
  59       a(2,npt)=-1./r
  60       a(5,npt)=-1./r
  61       npt=4+(j-1)*4
  62       a(1,npt)=3./r+1./(r/2.+rext)+ro/dt
  63       a(2,npt)=0.
  64       a(5,npt)=-1./r
  65       enddo
  66       i=1
  67       a(1,i)=2./r+ro/dt
  68       a(2,i)=-1./r
  69       a(5,i)=-1./r
  70       i=21
  71       a(1,i)=2./r+ro/dt
  72       a(2,i)=-1./r
  73       a(5,i)=-1./r
  74       i=24
  75       a(1,i)=2./r+1./(r/2.+rext)+ro/dt
  76       a(2,i)=-1./r
  77       a(5,i)=-1./r
  78       i=4
  79       a(1,i)=2./r+1./(r/2.+rext)+ro/dt
  80       a(2,i)=0.
  81       a(5,i)=-1./r
  82 c
  83 c  factorisation de la matrice
  84 c
  85       n=24
  86       kd=4
  87
  88       call dPBTRF( 'L' , N, KD, A, LDAB, INFO )
  89
  90        maille(1)=4
  91        maille(2)=6
  92
  93       ti=0.
  94
  95       do i=1,6
  96       tpar(i)=t(4*i)
  97       rpar(i)=r/2.
  98       enddo
  99 c
 100 c  initialisation de la temperature a t=0
 101 c
 102         CALL cpeRE(compo,CP_TEMPS, ti, npas, 'temperature', 24,
 103      &            t , info)
 104             IF( info.NE. CPOK  )GO TO 9000
 105 c
 106 c  initialisation de la temperature de bord a t=0.
 107 c
 108         CALL cpeRE(compo,CP_TEMPS, ti, npas, 'tparoi', 6,
 109      &            tpar , info)
 110             IF( info.NE. CPOK  )GO TO 9000
 111 c
 112 c  initialisation de la resistance thermique de bord a t=0.
 113 c
 114         CALL cpeRE(compo,CP_TEMPS, ti, npas, 'rparoi', 6,
 115      &            rpar , info)
 116             IF( info.NE. CPOK  )GO TO 9000
 117 c
 118 c  boucle temporelle infinie
 119 c
 120       do while( .TRUE. )
 121 c     do while( ti.lT.100. )
 122       tf=ti+dt
 123 c
 124 c  lecture de la puissance combustible entre ti et tf
 125 c
 126         CALL cplRE(compo,CP_TEMPS,ti, tf, npas, 'puissa', 24,
 127      &      nval, puiss , info)
 128             IF( info.NE. CPOK  )GO TO 9000
 129 c
 130 c  lecture de la temperature exterieure entre ti et tf
 131 c
 132         CALL cplRE(compo,CP_TEMPS,ti, tf, npas, 'text', 6,
 133      &      nval, text , info)
 134             IF( info.NE. CPOK  )GO TO 9000
 135 c
 136 c  lecture de la resistance exterieure entre ti et tf
 137 c
 138         CALL cplRE(compo,CP_TEMPS,ti, tf, npas, 'rext', 6,
 139      &      nval, rflu , info)
 140             IF( info.NE. CPOK  )GO TO 9000
 141 c
 142 c  calcul du second membre
 143 c
 144       do i=1,24
 145       b(i)=ro/dt*tn(i)
 146       enddo
 147
 148       do j=1,6
 149       npt=j*4
 150       b(npt)=b(npt)+text(j)/(r/2+rflu(j))
 151       enddo
 152
 153       do npt=1,24
 154       b(npt)=b(npt)+puiss(npt)
 155       enddo
 156 c
 157 c  resolution du systeme lineaire
 158 c
 159       call dPBTRs( 'L' , N,kd, nrhs,A,LDAB,b,ldb,INFO )
 160
 161       do i=1,24
 162       tn(i)=b(i)
 163       t(i)=b(i)
 164       enddo
 165
 166       do i=1,6
 167       tpar(i)=t(4*i)
 168       rpar(i)=r/2.
 169       enddo
 170       write(6,*)"SOLIDE: temps=",tf
 171       call flush(6)
 172 c
 173 c  ecriture de la temperature a t=tf
 174 c
 175         CALL cpeRE(compo,CP_TEMPS, tf, npas, 'temperature', 24,
 176      &            t , info)
 177             IF( info.NE. CPOK  )GO TO 9000
 178 c
 179 c  ecriture de la temperature de paroi a t=tf
 180 c
 181         CALL cpeRE(compo,CP_TEMPS, tf, npas, 'tparoi', 6,
 182      &            tpar , info)
 183             IF( info.NE. CPOK  )GO TO 9000
 184 c
 185 c  ecriture de la resistance de bord a t=tf
 186 c
 187         CALL cpeRE(compo,CP_TEMPS, tf, npas, 'rparoi', 6,
 188      &            rpar , info)
 189             IF( info.NE. CPOK  )GO TO 9000
 190
 191       ti=tf
 192
 193       enddo
 194 9000  continue
 195             CALL cpfin(compo,CP_ARRET, info)
 196       end
 197 c
 198       SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
 199 *
 200 *  -- LAPACK routine (version 2.0) --
 201 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 202 *     Courant Institute, Argonne National Lab, and Rice University
 203 *     March 31, 1993
 204 *
 205 *     .. Scalar Arguments ..
 206       CHARACTER          UPLO
 207       INTEGER            INFO, KD, LDAB, N
 208 *     ..
 209 *     .. Array Arguments ..
 210       DOUBLE PRECISION   AB( LDAB, * )
 211 *     ..
 212 *
 213 *  Purpose
 214 *  =======
 215 *
 216 *  DPBTRF computes the Cholesky factorization of a real symmetric
 217 *  positive definite band matrix A.
 218 *
 219 *  The factorization has the form
 220 *     A = U**T * U,  if UPLO = 'U', or
 221 *     A = L  * L**T,  if UPLO = 'L',
 222 *  where U is an upper triangular matrix and L is lower triangular.
 223 *
 224 *  Arguments
 225 *  =========
 226 *
 227 *  UPLO    (input) CHARACTER*1
 228 *          = 'U':  Upper triangle of A is stored;
 229 *          = 'L':  Lower triangle of A is stored.
 230 *
 231 *  N       (input) INTEGER
 232 *          The order of the matrix A.  N >= 0.
 233 *
 234 *  KD      (input) INTEGER
 235 *          The number of superdiagonals of the matrix A if UPLO = 'U',
 236 *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
 237 *
 238 *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
 239 *          On entry, the upper or lower triangle of the symmetric band
 240 *          matrix A, stored in the first KD+1 rows of the array.  The
 241 *          j-th column of A is stored in the j-th column of the array AB
 242 *          as follows:
 243 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
 244 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
 245 *
 246 *          On exit, if INFO = 0, the triangular factor U or L from the
 247 *          Cholesky factorization A = U**T*U or A = L*L**T of the band
 248 *          matrix A, in the same storage format as A.
 249 *
 250 *  LDAB    (input) INTEGER
 251 *          The leading dimension of the array AB.  LDAB >= KD+1.
 252 *
 253 *  INFO    (output) INTEGER
 254 *          = 0:  successful exit
 255 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 256 *          > 0:  if INFO = i, the leading minor of order i is not
 257 *                positive definite, and the factorization could not be
 258 *                completed.
 259 *
 260 *  Further Details
 261 *  ===============
 262 *
 263 *  The band storage scheme is illustrated by the following example, when
 264 *  N = 6, KD = 2, and UPLO = 'U':
 265 *
 266 *  On entry:                       On exit:
 267 *
 268 *      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
 269 *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
 270 *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
 271 *
 272 *  Similarly, if UPLO = 'L' the format of A is as follows:
 273 *
 274 *  On entry:                       On exit:
 275 *
 276 *     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
 277 *     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
 278 *     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
 279 *
 280 *  Array elements marked * are not used by the routine.
 281 *
 282 *  Contributed by
 283 *  Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
 284 *
 285 *  =====================================================================
 286 *
 287 *     .. Parameters ..
 288       DOUBLE PRECISION   ONE, ZERO
 289       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 290       INTEGER            NBMAX, LDWORK
 291       PARAMETER          ( NBMAX = 32, LDWORK = NBMAX+1 )
 292 *     ..
 293 *     .. Local Scalars ..
 294       INTEGER            I, I2, I3, IB, II, J, JJ, NB
 295 *     ..
 296 *     .. Local Arrays ..
 297       DOUBLE PRECISION   WORK( LDWORK, NBMAX )
 298 *     ..
 299 *     .. External Functions ..
 300       LOGICAL            LSAME
 301       INTEGER            ILAENV
 302       EXTERNAL           LSAME, ILAENV
 303 *     ..
 304 *     .. External Subroutines ..
 305       EXTERNAL           DGEMM, DPBTF2, DPOTF2, DSYRK, DTRSM, XERBLA
 306 *     ..
 307 *     .. Intrinsic Functions ..
 308       INTRINSIC          MIN
 309 *     ..
 310 *     .. Executable Statements ..
 311 *
 312 *     Test the input parameters.
 313 *
 314       INFO = 0
 315       IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND.
 316      $    ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
 317          INFO = -1
 318       ELSE IF( N.LT.0 ) THEN
 319          INFO = -2
 320       ELSE IF( KD.LT.0 ) THEN
 321          INFO = -3
 322       ELSE IF( LDAB.LT.KD+1 ) THEN
 323          INFO = -5
 324       END IF
 325       IF( INFO.NE.0 ) THEN
 326          CALL XERBLA( 'DPBTRF', -INFO )
 327          RETURN
 328       END IF
 329 *
 330 *     Quick return if possible
 331 *
 332       IF( N.EQ.0 )
 333      $   RETURN
 334 *
 335 *     Determine the block size for this environment
 336 *
 337       NB = ILAENV( 1, 'DPBTRF', UPLO, N, KD, -1, -1 )
 338 *
 339 *     The block size must not exceed the semi-bandwidth KD, and must not
 340 *     exceed the limit set by the size of the local array WORK.
 341 *
 342       NB = MIN( NB, NBMAX )
 343 *
 344       IF( NB.LE.1 .OR. NB.GT.KD ) THEN
 345 *
 346 *        Use unblocked code
 347 *
 348          CALL DPBTF2( UPLO, N, KD, AB, LDAB, INFO )
 349       ELSE
 350 *
 351 *        Use blocked code
 352 *
 353          IF( LSAME( UPLO, 'U' ) ) THEN
 354 *
 355 *           Compute the Cholesky factorization of a symmetric band
 356 *           matrix, given the upper triangle of the matrix in band
 357 *           storage.
 358 *
 359 *           Zero the upper triangle of the work array.
 360 *
 361             DO 20 J = 1, NB
 362                DO 10 I = 1, J - 1
 363                   WORK( I, J ) = ZERO
 364    10          CONTINUE
 365    20       CONTINUE
 366 *
 367 *           Process the band matrix one diagonal block at a time.
 368 *
 369             DO 70 I = 1, N, NB
 370                IB = MIN( NB, N-I+1 )
 371 *
 372 *              Factorize the diagonal block
 373 *
 374                CALL DPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II )
 375                IF( II.NE.0 ) THEN
 376                   INFO = I + II - 1
 377                   GO TO 150
 378                END IF
 379                IF( I+IB.LE.N ) THEN
 380 *
 381 *                 Update the relevant part of the trailing submatrix.
 382 *                 If A11 denotes the diagonal block which has just been
 383 *                 factorized, then we need to update the remaining
 384 *                 blocks in the diagram:
 385 *
 386 *                    A11   A12   A13
 387 *                          A22   A23
 388 *                                A33
 389 *
 390 *                 The numbers of rows and columns in the partitioning
 391 *                 are IB, I2, I3 respectively. The blocks A12, A22 and
 392 *                 A23 are empty if IB = KD. The upper triangle of A13
 393 *                 lies outside the band.
 394 *
 395                   I2 = MIN( KD-IB, N-I-IB+1 )
 396                   I3 = MIN( IB, N-I-KD+1 )
 397 *
 398                   IF( I2.GT.0 ) THEN
 399 *
 400 *                    Update A12
 401 *
 402                      CALL DTRSM( 'Left', 'Upper', 'Transpose',
 403      $                           'Non-unit', IB, I2, ONE, AB( KD+1, I ),
 404      $                           LDAB-1, AB( KD+1-IB, I+IB ), LDAB-1 )
 405 *
 406 *                    Update A22
 407 *
 408                      CALL DSYRK( 'Upper', 'Transpose', I2, IB, -ONE,
 409      $                           AB( KD+1-IB, I+IB ), LDAB-1, ONE,
 410      $                           AB( KD+1, I+IB ), LDAB-1 )
 411                   END IF
 412 *
 413                   IF( I3.GT.0 ) THEN
 414 *
 415 *                    Copy the lower triangle of A13 into the work array.
 416 *
 417                      DO 40 JJ = 1, I3
 418                         DO 30 II = JJ, IB
 419                            WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 )
 420    30                   CONTINUE
 421    40                CONTINUE
 422 *
 423 *                    Update A13 (in the work array).
 424 *
 425                      CALL DTRSM( 'Left', 'Upper', 'Transpose',
 426      $                           'Non-unit', IB, I3, ONE, AB( KD+1, I ),
 427      $                           LDAB-1, WORK, LDWORK )
 428 *
 429 *                    Update A23
 430 *
 431                      IF( I2.GT.0 )
 432      $                  CALL DGEMM( 'Transpose', 'No Transpose', I2, I3,
 433      $                              IB, -ONE, AB( KD+1-IB, I+IB ),
 434      $                              LDAB-1, WORK, LDWORK, ONE,
 435      $                              AB( 1+IB, I+KD ), LDAB-1 )
 436 *
 437 *                    Update A33
 438 *
 439                      CALL DSYRK( 'Upper', 'Transpose', I3, IB, -ONE,
 440      $                           WORK, LDWORK, ONE, AB( KD+1, I+KD ),
 441      $                           LDAB-1 )
 442 *
 443 *                    Copy the lower triangle of A13 back into place.
 444 *
 445                      DO 60 JJ = 1, I3
 446                         DO 50 II = JJ, IB
 447                            AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ )
 448    50                   CONTINUE
 449    60                CONTINUE
 450                   END IF
 451                END IF
 452    70       CONTINUE
 453          ELSE
 454 *
 455 *           Compute the Cholesky factorization of a symmetric band
 456 *           matrix, given the lower triangle of the matrix in band
 457 *           storage.
 458 *
 459 *           Zero the lower triangle of the work array.
 460 *
 461             DO 90 J = 1, NB
 462                DO 80 I = J + 1, NB
 463                   WORK( I, J ) = ZERO
 464    80          CONTINUE
 465    90       CONTINUE
 466 *
 467 *           Process the band matrix one diagonal block at a time.
 468 *
 469             DO 140 I = 1, N, NB
 470                IB = MIN( NB, N-I+1 )
 471 *
 472 *              Factorize the diagonal block
 473 *
 474                CALL DPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II )
 475                IF( II.NE.0 ) THEN
 476                   INFO = I + II - 1
 477                   GO TO 150
 478                END IF
 479                IF( I+IB.LE.N ) THEN
 480 *
 481 *                 Update the relevant part of the trailing submatrix.
 482 *                 If A11 denotes the diagonal block which has just been
 483 *                 factorized, then we need to update the remaining
 484 *                 blocks in the diagram:
 485 *
 486 *                    A11
 487 *                    A21   A22
 488 *                    A31   A32   A33
 489 *
 490 *                 The numbers of rows and columns in the partitioning
 491 *                 are IB, I2, I3 respectively. The blocks A21, A22 and
 492 *                 A32 are empty if IB = KD. The lower triangle of A31
 493 *                 lies outside the band.
 494 *
 495                   I2 = MIN( KD-IB, N-I-IB+1 )
 496                   I3 = MIN( IB, N-I-KD+1 )
 497 *
 498                   IF( I2.GT.0 ) THEN
 499 *
 500 *                    Update A21
 501 *
 502                      CALL DTRSM( 'Right', 'Lower', 'Transpose',
 503      $                           'Non-unit', I2, IB, ONE, AB( 1, I ),
 504      $                           LDAB-1, AB( 1+IB, I ), LDAB-1 )
 505 *
 506 *                    Update A22
 507 *
 508                      CALL DSYRK( 'Lower', 'No Transpose', I2, IB, -ONE,
 509      $                           AB( 1+IB, I ), LDAB-1, ONE,
 510      $                           AB( 1, I+IB ), LDAB-1 )
 511                   END IF
 512 *
 513                   IF( I3.GT.0 ) THEN
 514 *
 515 *                    Copy the upper triangle of A31 into the work array.
 516 *
 517                      DO 110 JJ = 1, IB
 518                         DO 100 II = 1, MIN( JJ, I3 )
 519                            WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 )
 520   100                   CONTINUE
 521   110                CONTINUE
 522 *
 523 *                    Update A31 (in the work array).
 524 *
 525                      CALL DTRSM( 'Right', 'Lower', 'Transpose',
 526      $                           'Non-unit', I3, IB, ONE, AB( 1, I ),
 527      $                           LDAB-1, WORK, LDWORK )
 528 *
 529 *                    Update A32
 530 *
 531                      IF( I2.GT.0 )
 532      $                  CALL DGEMM( 'No transpose', 'Transpose', I3, I2,
 533      $                              IB, -ONE, WORK, LDWORK,
 534      $                              AB( 1+IB, I ), LDAB-1, ONE,
 535      $                              AB( 1+KD-IB, I+IB ), LDAB-1 )
 536 *
 537 *                    Update A33
 538 *
 539                      CALL DSYRK( 'Lower', 'No Transpose', I3, IB, -ONE,
 540      $                           WORK, LDWORK, ONE, AB( 1, I+KD ),
 541      $                           LDAB-1 )
 542 *
 543 *                    Copy the upper triangle of A31 back into place.
 544 *
 545                      DO 130 JJ = 1, IB
 546                         DO 120 II = 1, MIN( JJ, I3 )
 547                            AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ )
 548   120                   CONTINUE
 549   130                CONTINUE
 550                   END IF
 551                END IF
 552   140       CONTINUE
 553          END IF
 554       END IF
 555       RETURN
 556 *
 557   150 CONTINUE
 558       RETURN
 559 *
 560 *     End of DPBTRF
 561 *
 562       END
 563       SUBROUTINE DPBTRS( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO )
 564 *
 565 *  -- LAPACK routine (version 2.0) --
 566 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 567 *     Courant Institute, Argonne National Lab, and Rice University
 568 *     September 30, 1994
 569 *
 570 *     .. Scalar Arguments ..
 571       CHARACTER          UPLO
 572       INTEGER            INFO, KD, LDAB, LDB, N, NRHS
 573 *     ..
 574 *     .. Array Arguments ..
 575       DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * )
 576 *     ..
 577 *
 578 *  Purpose
 579 *  =======
 580 *
 581 *  DPBTRS solves a system of linear equations A*X = B with a symmetric
 582 *  positive definite band matrix A using the Cholesky factorization
 583 *  A = U**T*U or A = L*L**T computed by DPBTRF.
 584 *
 585 *  Arguments
 586 *  =========
 587 *
 588 *  UPLO    (input) CHARACTER*1
 589 *          = 'U':  Upper triangular factor stored in AB;
 590 *          = 'L':  Lower triangular factor stored in AB.
 591 *
 592 *  N       (input) INTEGER
 593 *          The order of the matrix A.  N >= 0.
 594 *
 595 *  KD      (input) INTEGER
 596 *          The number of superdiagonals of the matrix A if UPLO = 'U',
 597 *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
 598 *
 599 *  NRHS    (input) INTEGER
 600 *          The number of right hand sides, i.e., the number of columns
 601 *          of the matrix B.  NRHS >= 0.
 602 *
 603 *  AB      (input) DOUBLE PRECISION array, dimension (LDAB,N)
 604 *          The triangular factor U or L from the Cholesky factorization
 605 *          A = U**T*U or A = L*L**T of the band matrix A, stored in the
 606 *          first KD+1 rows of the array.  The j-th column of U or L is
 607 *          stored in the j-th column of the array AB as follows:
 608 *          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
 609 *          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).
 610 *
 611 *  LDAB    (input) INTEGER
 612 *          The leading dimension of the array AB.  LDAB >= KD+1.
 613 *
 614 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 615 *          On entry, the right hand side matrix B.
 616 *          On exit, the solution matrix X.
 617 *
 618 *  LDB     (input) INTEGER
 619 *          The leading dimension of the array B.  LDB >= max(1,N).
 620 *
 621 *  INFO    (output) INTEGER
 622 *          = 0:  successful exit
 623 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 624 *
 625 *  =====================================================================
 626 *
 627 *     .. Local Scalars ..
 628       LOGICAL            UPPER
 629       INTEGER            J
 630 *     ..
 631 *     .. External Functions ..
 632       LOGICAL            LSAME
 633       EXTERNAL           LSAME
 634 *     ..
 635 *     .. External Subroutines ..
 636       EXTERNAL           DTBSV, XERBLA
 637 *     ..
 638 *     .. Intrinsic Functions ..
 639       INTRINSIC          MAX
 640 *     ..
 641 *     .. Executable Statements ..
 642 *
 643 *     Test the input parameters.
 644 *
 645       INFO = 0
 646       UPPER = LSAME( UPLO, 'U' )
 647       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 648          INFO = -1
 649       ELSE IF( N.LT.0 ) THEN
 650          INFO = -2
 651       ELSE IF( KD.LT.0 ) THEN
 652          INFO = -3
 653       ELSE IF( NRHS.LT.0 ) THEN
 654          INFO = -4
 655       ELSE IF( LDAB.LT.KD+1 ) THEN
 656          INFO = -6
 657       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 658          INFO = -8
 659       END IF
 660       IF( INFO.NE.0 ) THEN
 661          CALL XERBLA( 'DPBTRS', -INFO )
 662          RETURN
 663       END IF
 664 *
 665 *     Quick return if possible
 666 *
 667       IF( N.EQ.0 .OR. NRHS.EQ.0 )
 668      $   RETURN
 669 *
 670       IF( UPPER ) THEN
 671 *
 672 *        Solve A*X = B where A = U'*U.
 673 *
 674          DO 10 J = 1, NRHS
 675 *
 676 *           Solve U'*X = B, overwriting B with X.
 677 *
 678             CALL DTBSV( 'Upper', 'Transpose', 'Non-unit', N, KD, AB,
 679      $                  LDAB, B( 1, J ), 1 )
 680 *
 681 *           Solve U*X = B, overwriting B with X.
 682 *
 683             CALL DTBSV( 'Upper', 'No transpose', 'Non-unit', N, KD, AB,
 684      $                  LDAB, B( 1, J ), 1 )
 685    10    CONTINUE
 686       ELSE
 687 *
 688 *        Solve A*X = B where A = L*L'.
 689 *
 690          DO 20 J = 1, NRHS
 691 *
 692 *           Solve L*X = B, overwriting B with X.
 693 *
 694             CALL DTBSV( 'Lower', 'No transpose', 'Non-unit', N, KD, AB,
 695      $                  LDAB, B( 1, J ), 1 )
 696 *
 697 *           Solve L'*X = B, overwriting B with X.
 698 *
 699             CALL DTBSV( 'Lower', 'Transpose', 'Non-unit', N, KD, AB,
 700      $                  LDAB, B( 1, J ), 1 )
 701    20    CONTINUE
 702       END IF
 703 *
 704       RETURN
 705 *
 706 *     End of DPBTRS
 707 *
 708       END
 709       INTEGER          FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,
 710      $                 N4 )
 711 *
 712 *  -- LAPACK auxiliary routine (version 2.0) --
 713 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 714 *     Courant Institute, Argonne National Lab, and Rice University
 715 *     September 30, 1994
 716 *
 717 *     .. Scalar Arguments ..
 718       CHARACTER*( * )    NAME, OPTS
 719       INTEGER            ISPEC, N1, N2, N3, N4
 720 *     ..
 721 *
 722 *  Purpose
 723 *  =======
 724 *
 725 *  ILAENV is called from the LAPACK routines to choose problem-dependent
 726 *  parameters for the local environment.  See ISPEC for a description of
 727 *  the parameters.
 728 *
 729 *  This version provides a set of parameters which should give good,
 730 *  but not optimal, performance on many of the currently available
 731 *  computers.  Users are encouraged to modify this subroutine to set
 732 *  the tuning parameters for their particular machine using the option
 733 *  and problem size information in the arguments.
 734 *
 735 *  This routine will not function correctly if it is converted to all
 736 *  lower case.  Converting it to all upper case is allowed.
 737 *
 738 *  Arguments
 739 *  =========
 740 *
 741 *  ISPEC   (input) INTEGER
 742 *          Specifies the parameter to be returned as the value of
 743 *          ILAENV.
 744 *          = 1: the optimal blocksize; if this value is 1, an unblocked
 745 *               algorithm will give the best performance.
 746 *          = 2: the minimum block size for which the block routine
 747 *               should be used; if the usable block size is less than
 748 *               this value, an unblocked routine should be used.
 749 *          = 3: the crossover point (in a block routine, for N less
 750 *               than this value, an unblocked routine should be used)
 751 *          = 4: the number of shifts, used in the nonsymmetric
 752 *               eigenvalue routines
 753 *          = 5: the minimum column dimension for blocking to be used;
 754 *               rectangular blocks must have dimension at least k by m,
 755 *               where k is given by ILAENV(2,...) and m by ILAENV(5,...)
 756 *          = 6: the crossover point for the SVD (when reducing an m by n
 757 *               matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
 758 *               this value, a QR factorization is used first to reduce
 759 *               the matrix to a triangular form.)
 760 *          = 7: the number of processors
 761 *          = 8: the crossover point for the multishift QR and QZ methods
 762 *               for nonsymmetric eigenvalue problems.
 763 *
 764 *  NAME    (input) CHARACTER*(*)
 765 *          The name of the calling subroutine, in either upper case or
 766 *          lower case.
 767 *
 768 *  OPTS    (input) CHARACTER*(*)
 769 *          The character options to the subroutine NAME, concatenated
 770 *          into a single character string.  For example, UPLO = 'U',
 771 *          TRANS = 'T', and DIAG = 'N' for a triangular routine would
 772 *          be specified as OPTS = 'UTN'.
 773 *
 774 *  N1      (input) INTEGER
 775 *  N2      (input) INTEGER
 776 *  N3      (input) INTEGER
 777 *  N4      (input) INTEGER
 778 *          Problem dimensions for the subroutine NAME; these may not all
 779 *          be required.
 780 *
 781 * (ILAENV) (output) INTEGER
 782 *          >= 0: the value of the parameter specified by ISPEC
 783 *          < 0:  if ILAENV = -k, the k-th argument had an illegal value.
 784 *
 785 *  Further Details
 786 *  ===============
 787 *
 788 *  The following conventions have been used when calling ILAENV from the
 789 *  LAPACK routines:
 790 *  1)  OPTS is a concatenation of all of the character options to
 791 *      subroutine NAME, in the same order that they appear in the
 792 *      argument list for NAME, even if they are not used in determining
 793 *      the value of the parameter specified by ISPEC.
 794 *  2)  The problem dimensions N1, N2, N3, N4 are specified in the order
 795 *      that they appear in the argument list for NAME.  N1 is used
 796 *      first, N2 second, and so on, and unused problem dimensions are
 797 *      passed a value of -1.
 798 *  3)  The parameter value returned by ILAENV is checked for validity in
 799 *      the calling subroutine.  For example, ILAENV is used to retrieve
 800 *      the optimal blocksize for STRTRI as follows:
 801 *
 802 *      NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
 803 *      IF( NB.LE.1 ) NB = MAX( 1, N )
 804 *
 805 *  =====================================================================
 806 *
 807 *     .. Local Scalars ..
 808       LOGICAL            CNAME, SNAME
 809       CHARACTER*1        C1
 810       CHARACTER*2        C2, C4
 811       CHARACTER*3        C3
 812       CHARACTER*6        SUBNAM
 813       INTEGER            I, IC, IZ, NB, NBMIN, NX
 814 *     ..
 815 *     .. Intrinsic Functions ..
 816       INTRINSIC          CHAR, ICHAR, INT, MIN, REAL
 817 *     ..
 818 *     .. Executable Statements ..
 819 *
 820       GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC
 821 *
 822 *     Invalid value for ISPEC
 823 *
 824       ILAENV = -1
 825       RETURN
 826 *
 827   100 CONTINUE
 828 *
 829 *     Convert NAME to upper case if the first character is lower case.
 830 *
 831       ILAENV = 1
 832       SUBNAM = NAME
 833       IC = ICHAR( SUBNAM( 1:1 ) )
 834       IZ = ICHAR( 'Z' )
 835       IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN
 836 *
 837 *        ASCII character set
 838 *
 839          IF( IC.GE.97 .AND. IC.LE.122 ) THEN
 840             SUBNAM( 1:1 ) = CHAR( IC-32 )
 841             DO 10 I = 2, 6
 842                IC = ICHAR( SUBNAM( I:I ) )
 843                IF( IC.GE.97 .AND. IC.LE.122 )
 844      $            SUBNAM( I:I ) = CHAR( IC-32 )
 845    10       CONTINUE
 846          END IF
 847 *
 848       ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN
 849 *
 850 *        EBCDIC character set
 851 *
 852          IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
 853      $       ( IC.GE.145 .AND. IC.LE.153 ) .OR.
 854      $       ( IC.GE.162 .AND. IC.LE.169 ) ) THEN
 855             SUBNAM( 1:1 ) = CHAR( IC+64 )
 856             DO 20 I = 2, 6
 857                IC = ICHAR( SUBNAM( I:I ) )
 858                IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
 859      $             ( IC.GE.145 .AND. IC.LE.153 ) .OR.
 860      $             ( IC.GE.162 .AND. IC.LE.169 ) )
 861      $            SUBNAM( I:I ) = CHAR( IC+64 )
 862    20       CONTINUE
 863          END IF
 864 *
 865       ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN
 866 *
 867 *        Prime machines:  ASCII+128
 868 *
 869          IF( IC.GE.225 .AND. IC.LE.250 ) THEN
 870             SUBNAM( 1:1 ) = CHAR( IC-32 )
 871             DO 30 I = 2, 6
 872                IC = ICHAR( SUBNAM( I:I ) )
 873                IF( IC.GE.225 .AND. IC.LE.250 )
 874      $            SUBNAM( I:I ) = CHAR( IC-32 )
 875    30       CONTINUE
 876          END IF
 877       END IF
 878 *
 879       C1 = SUBNAM( 1:1 )
 880       SNAME = C1.EQ.'S' .OR. C1.EQ.'D'
 881       CNAME = C1.EQ.'C' .OR. C1.EQ.'Z'
 882       IF( .NOT.( CNAME .OR. SNAME ) )
 883      $   RETURN
 884       C2 = SUBNAM( 2:3 )
 885       C3 = SUBNAM( 4:6 )
 886       C4 = C3( 2:3 )
 887 *
 888       GO TO ( 110, 200, 300 ) ISPEC
 889 *
 890   110 CONTINUE
 891 *
 892 *     ISPEC = 1:  block size
 893 *
 894 *     In these examples, separate code is provided for setting NB for
 895 *     real and complex.  We assume that NB will take the same value in
 896 *     single or double precision.
 897 *
 898       NB = 1
 899 *
 900       IF( C2.EQ.'GE' ) THEN
 901          IF( C3.EQ.'TRF' ) THEN
 902             IF( SNAME ) THEN
 903                NB = 64
 904             ELSE
 905                NB = 64
 906             END IF
 907          ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
 908      $            C3.EQ.'QLF' ) THEN
 909             IF( SNAME ) THEN
 910                NB = 32
 911             ELSE
 912                NB = 32
 913             END IF
 914          ELSE IF( C3.EQ.'HRD' ) THEN
 915             IF( SNAME ) THEN
 916                NB = 32
 917             ELSE
 918                NB = 32
 919             END IF
 920          ELSE IF( C3.EQ.'BRD' ) THEN
 921             IF( SNAME ) THEN
 922                NB = 32
 923             ELSE
 924                NB = 32
 925             END IF
 926          ELSE IF( C3.EQ.'TRI' ) THEN
 927             IF( SNAME ) THEN
 928                NB = 64
 929             ELSE
 930                NB = 64
 931             END IF
 932          END IF
 933       ELSE IF( C2.EQ.'PO' ) THEN
 934          IF( C3.EQ.'TRF' ) THEN
 935             IF( SNAME ) THEN
 936                NB = 64
 937             ELSE
 938                NB = 64
 939             END IF
 940          END IF
 941       ELSE IF( C2.EQ.'SY' ) THEN
 942          IF( C3.EQ.'TRF' ) THEN
 943             IF( SNAME ) THEN
 944                NB = 64
 945             ELSE
 946                NB = 64
 947             END IF
 948          ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
 949             NB = 1
 950          ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN
 951             NB = 64
 952          END IF
 953       ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
 954          IF( C3.EQ.'TRF' ) THEN
 955             NB = 64
 956          ELSE IF( C3.EQ.'TRD' ) THEN
 957             NB = 1
 958          ELSE IF( C3.EQ.'GST' ) THEN
 959             NB = 64
 960          END IF
 961       ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
 962          IF( C3( 1:1 ).EQ.'G' ) THEN
 963             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
 964      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
 965      $          C4.EQ.'BR' ) THEN
 966                NB = 32
 967             END IF
 968          ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
 969             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
 970      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
 971      $          C4.EQ.'BR' ) THEN
 972                NB = 32
 973             END IF
 974          END IF
 975       ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
 976          IF( C3( 1:1 ).EQ.'G' ) THEN
 977             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
 978      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
 979      $          C4.EQ.'BR' ) THEN
 980                NB = 32
 981             END IF
 982          ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
 983             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
 984      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
 985      $          C4.EQ.'BR' ) THEN
 986                NB = 32
 987             END IF
 988          END IF
 989       ELSE IF( C2.EQ.'GB' ) THEN
 990          IF( C3.EQ.'TRF' ) THEN
 991             IF( SNAME ) THEN
 992                IF( N4.LE.64 ) THEN
 993                   NB = 1
 994                ELSE
 995                   NB = 32
 996                END IF
 997             ELSE
 998                IF( N4.LE.64 ) THEN
 999                   NB = 1
1000                ELSE
1001                   NB = 32
1002                END IF
1003             END IF
1004          END IF
1005       ELSE IF( C2.EQ.'PB' ) THEN
1006          IF( C3.EQ.'TRF' ) THEN
1007             IF( SNAME ) THEN
1008                IF( N2.LE.64 ) THEN
1009                   NB = 1
1010                ELSE
1011                   NB = 32
1012                END IF
1013             ELSE
1014                IF( N2.LE.64 ) THEN
1015                   NB = 1
1016                ELSE
1017                   NB = 32
1018                END IF
1019             END IF
1020          END IF
1021       ELSE IF( C2.EQ.'TR' ) THEN
1022          IF( C3.EQ.'TRI' ) THEN
1023             IF( SNAME ) THEN
1024                NB = 64
1025             ELSE
1026                NB = 64
1027             END IF
1028          END IF
1029       ELSE IF( C2.EQ.'LA' ) THEN
1030          IF( C3.EQ.'UUM' ) THEN
1031             IF( SNAME ) THEN
1032                NB = 64
1033             ELSE
1034                NB = 64
1035             END IF
1036          END IF
1037       ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
1038          IF( C3.EQ.'EBZ' ) THEN
1039             NB = 1
1040          END IF
1041       END IF
1042       ILAENV = NB
1043       RETURN
1044 *
1045   200 CONTINUE
1046 *
1047 *     ISPEC = 2:  minimum block size
1048 *
1049       NBMIN = 2
1050       IF( C2.EQ.'GE' ) THEN
1051          IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
1052      $       C3.EQ.'QLF' ) THEN
1053             IF( SNAME ) THEN
1054                NBMIN = 2
1055             ELSE
1056                NBMIN = 2
1057             END IF
1058          ELSE IF( C3.EQ.'HRD' ) THEN
1059             IF( SNAME ) THEN
1060                NBMIN = 2
1061             ELSE
1062                NBMIN = 2
1063             END IF
1064          ELSE IF( C3.EQ.'BRD' ) THEN
1065             IF( SNAME ) THEN
1066                NBMIN = 2
1067             ELSE
1068                NBMIN = 2
1069             END IF
1070          ELSE IF( C3.EQ.'TRI' ) THEN
1071             IF( SNAME ) THEN
1072                NBMIN = 2
1073             ELSE
1074                NBMIN = 2
1075             END IF
1076          END IF
1077       ELSE IF( C2.EQ.'SY' ) THEN
1078          IF( C3.EQ.'TRF' ) THEN
1079             IF( SNAME ) THEN
1080                NBMIN = 8
1081             ELSE
1082                NBMIN = 8
1083             END IF
1084          ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
1085             NBMIN = 2
1086          END IF
1087       ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
1088          IF( C3.EQ.'TRD' ) THEN
1089             NBMIN = 2
1090          END IF
1091       ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
1092          IF( C3( 1:1 ).EQ.'G' ) THEN
1093             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1094      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1095      $          C4.EQ.'BR' ) THEN
1096                NBMIN = 2
1097             END IF
1098          ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
1099             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1100      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1101      $          C4.EQ.'BR' ) THEN
1102                NBMIN = 2
1103             END IF
1104          END IF
1105       ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
1106          IF( C3( 1:1 ).EQ.'G' ) THEN
1107             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1108      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1109      $          C4.EQ.'BR' ) THEN
1110                NBMIN = 2
1111             END IF
1112          ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
1113             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1114      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1115      $          C4.EQ.'BR' ) THEN
1116                NBMIN = 2
1117             END IF
1118          END IF
1119       END IF
1120       ILAENV = NBMIN
1121       RETURN
1122 *
1123   300 CONTINUE
1124 *
1125 *     ISPEC = 3:  crossover point
1126 *
1127       NX = 0
1128       IF( C2.EQ.'GE' ) THEN
1129          IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
1130      $       C3.EQ.'QLF' ) THEN
1131             IF( SNAME ) THEN
1132                NX = 128
1133             ELSE
1134                NX = 128
1135             END IF
1136          ELSE IF( C3.EQ.'HRD' ) THEN
1137             IF( SNAME ) THEN
1138                NX = 128
1139             ELSE
1140                NX = 128
1141             END IF
1142          ELSE IF( C3.EQ.'BRD' ) THEN
1143             IF( SNAME ) THEN
1144                NX = 128
1145             ELSE
1146                NX = 128
1147             END IF
1148          END IF
1149       ELSE IF( C2.EQ.'SY' ) THEN
1150          IF( SNAME .AND. C3.EQ.'TRD' ) THEN
1151             NX = 1
1152          END IF
1153       ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
1154          IF( C3.EQ.'TRD' ) THEN
1155             NX = 1
1156          END IF
1157       ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
1158          IF( C3( 1:1 ).EQ.'G' ) THEN
1159             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1160      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1161      $          C4.EQ.'BR' ) THEN
1162                NX = 128
1163             END IF
1164          END IF
1165       ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
1166          IF( C3( 1:1 ).EQ.'G' ) THEN
1167             IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1168      $          C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1169      $          C4.EQ.'BR' ) THEN
1170                NX = 128
1171             END IF
1172          END IF
1173       END IF
1174       ILAENV = NX
1175       RETURN
1176 *
1177   400 CONTINUE
1178 *
1179 *     ISPEC = 4:  number of shifts (used by xHSEQR)
1180 *
1181       ILAENV = 6
1182       RETURN
1183 *
1184   500 CONTINUE
1185 *
1186 *     ISPEC = 5:  minimum column dimension (not used)
1187 *
1188       ILAENV = 2
1189       RETURN
1190 *
1191   600 CONTINUE
1192 *
1193 *     ISPEC = 6:  crossover point for SVD (used by xGELSS and xGESVD)
1194 *
1195       ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 )
1196       RETURN
1197 *
1198   700 CONTINUE
1199 *
1200 *     ISPEC = 7:  number of processors (not used)
1201 *
1202       ILAENV = 1
1203       RETURN
1204 *
1205   800 CONTINUE
1206 *
1207 *     ISPEC = 8:  crossover point for multishift (used by xHSEQR)
1208 *
1209       ILAENV = 50
1210       RETURN
1211 *
1212 *     End of ILAENV
1213 *
1214       END
1215       SUBROUTINE DPBTF2( UPLO, N, KD, AB, LDAB, INFO )
1216 *
1217 *  -- LAPACK routine (version 2.0) --
1218 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1219 *     Courant Institute, Argonne National Lab, and Rice University
1220 *     February 29, 1992
1221 *
1222 *     .. Scalar Arguments ..
1223       CHARACTER          UPLO
1224       INTEGER            INFO, KD, LDAB, N
1225 *     ..
1226 *     .. Array Arguments ..
1227       DOUBLE PRECISION   AB( LDAB, * )
1228 *     ..
1229 *
1230 *  Purpose
1231 *  =======
1232 *
1233 *  DPBTF2 computes the Cholesky factorization of a real symmetric
1234 *  positive definite band matrix A.
1235 *
1236 *  The factorization has the form
1237 *     A = U' * U ,  if UPLO = 'U', or
1238 *     A = L  * L',  if UPLO = 'L',
1239 *  where U is an upper triangular matrix, U' is the transpose of U, and
1240 *  L is lower triangular.
1241 *
1242 *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
1243 *
1244 *  Arguments
1245 *  =========
1246 *
1247 *  UPLO    (input) CHARACTER*1
1248 *          Specifies whether the upper or lower triangular part of the
1249 *          symmetric matrix A is stored:
1250 *          = 'U':  Upper triangular
1251 *          = 'L':  Lower triangular
1252 *
1253 *  N       (input) INTEGER
1254 *          The order of the matrix A.  N >= 0.
1255 *
1256 *  KD      (input) INTEGER
1257 *          The number of super-diagonals of the matrix A if UPLO = 'U',
1258 *          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
1259 *
1260 *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
1261 *          On entry, the upper or lower triangle of the symmetric band
1262 *          matrix A, stored in the first KD+1 rows of the array.  The
1263 *          j-th column of A is stored in the j-th column of the array AB
1264 *          as follows:
1265 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
1266 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
1267 *
1268 *          On exit, if INFO = 0, the triangular factor U or L from the
1269 *          Cholesky factorization A = U'*U or A = L*L' of the band
1270 *          matrix A, in the same storage format as A.
1271 *
1272 *  LDAB    (input) INTEGER
1273 *          The leading dimension of the array AB.  LDAB >= KD+1.
1274 *
1275 *  INFO    (output) INTEGER
1276 *          = 0: successful exit
1277 *          < 0: if INFO = -k, the k-th argument had an illegal value
1278 *          > 0: if INFO = k, the leading minor of order k is not
1279 *               positive definite, and the factorization could not be
1280 *               completed.
1281 *
1282 *  Further Details
1283 *  ===============
1284 *
1285 *  The band storage scheme is illustrated by the following example, when
1286 *  N = 6, KD = 2, and UPLO = 'U':
1287 *
1288 *  On entry:                       On exit:
1289 *
1290 *      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
1291 *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
1292 *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
1293 *
1294 *  Similarly, if UPLO = 'L' the format of A is as follows:
1295 *
1296 *  On entry:                       On exit:
1297 *
1298 *     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
1299 *     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
1300 *     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
1301 *
1302 *  Array elements marked * are not used by the routine.
1303 *
1304 *  =====================================================================
1305 *
1306 *     .. Parameters ..
1307       DOUBLE PRECISION   ONE, ZERO
1308       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
1309 *     ..
1310 *     .. Local Scalars ..
1311       LOGICAL            UPPER
1312       INTEGER            J, KLD, KN
1313       DOUBLE PRECISION   AJJ
1314 *     ..
1315 *     .. External Functions ..
1316       LOGICAL            LSAME
1317       EXTERNAL           LSAME
1318 *     ..
1319 *     .. External Subroutines ..
1320       EXTERNAL           DSCAL, DSYR, XERBLA
1321 *     ..
1322 *     .. Intrinsic Functions ..
1323       INTRINSIC          MAX, MIN, SQRT
1324 *     ..
1325 *     .. Executable Statements ..
1326 *
1327 *     Test the input parameters.
1328 *
1329       INFO = 0
1330       UPPER = LSAME( UPLO, 'U' )
1331       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
1332          INFO = -1
1333       ELSE IF( N.LT.0 ) THEN
1334          INFO = -2
1335       ELSE IF( KD.LT.0 ) THEN
1336          INFO = -3
1337       ELSE IF( LDAB.LT.KD+1 ) THEN
1338          INFO = -5
1339       END IF
1340       IF( INFO.NE.0 ) THEN
1341          CALL XERBLA( 'DPBTF2', -INFO )
1342          RETURN
1343       END IF
1344 *
1345 *     Quick return if possible
1346 *
1347       IF( N.EQ.0 )
1348      $   RETURN
1349 *
1350       KLD = MAX( 1, LDAB-1 )
1351 *
1352       IF( UPPER ) THEN
1353 *
1354 *        Compute the Cholesky factorization A = U'*U.
1355 *
1356          DO 10 J = 1, N
1357 *
1358 *           Compute U(J,J) and test for non-positive-definiteness.
1359 *
1360             AJJ = AB( KD+1, J )
1361             IF( AJJ.LE.ZERO )
1362      $         GO TO 30
1363             AJJ = SQRT( AJJ )
1364             AB( KD+1, J ) = AJJ
1365 *
1366 *           Compute elements J+1:J+KN of row J and update the
1367 *           trailing submatrix within the band.
1368 *
1369             KN = MIN( KD, N-J )
1370             IF( KN.GT.0 ) THEN
1371                CALL DSCAL( KN, ONE / AJJ, AB( KD, J+1 ), KLD )
1372                CALL DSYR( 'Upper', KN, -ONE, AB( KD, J+1 ), KLD,
1373      $                    AB( KD+1, J+1 ), KLD )
1374             END IF
1375    10    CONTINUE
1376       ELSE
1377 *
1378 *        Compute the Cholesky factorization A = L*L'.
1379 *
1380          DO 20 J = 1, N
1381 *
1382 *           Compute L(J,J) and test for non-positive-definiteness.
1383 *
1384             AJJ = AB( 1, J )
1385             IF( AJJ.LE.ZERO )
1386      $         GO TO 30
1387             AJJ = SQRT( AJJ )
1388             AB( 1, J ) = AJJ
1389 *
1390 *           Compute elements J+1:J+KN of column J and update the
1391 *           trailing submatrix within the band.
1392 *
1393             KN = MIN( KD, N-J )
1394             IF( KN.GT.0 ) THEN
1395                CALL DSCAL( KN, ONE / AJJ, AB( 2, J ), 1 )
1396                CALL DSYR( 'Lower', KN, -ONE, AB( 2, J ), 1,
1397      $                    AB( 1, J+1 ), KLD )
1398             END IF
1399    20    CONTINUE
1400       END IF
1401       RETURN
1402 *
1403    30 CONTINUE
1404       INFO = J
1405       RETURN
1406 *
1407 *     End of DPBTF2
1408 *
1409       END
1410       LOGICAL          FUNCTION LSAME( CA, CB )
1411 *
1412 *  -- LAPACK auxiliary routine (version 2.0) --
1413 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1414 *     Courant Institute, Argonne National Lab, and Rice University
1415 *     September 30, 1994
1416 *
1417 *     .. Scalar Arguments ..
1418       CHARACTER          CA, CB
1419 *     ..
1420 *
1421 *  Purpose
1422 *  =======
1423 *
1424 *  LSAME returns .TRUE. if CA is the same letter as CB regardless of
1425 *  case.
1426 *
1427 *  Arguments
1428 *  =========
1429 *
1430 *  CA      (input) CHARACTER*1
1431 *  CB      (input) CHARACTER*1
1432 *          CA and CB specify the single characters to be compared.
1433 *
1434 * =====================================================================
1435 *
1436 *     .. Intrinsic Functions ..
1437       INTRINSIC          ICHAR
1438 *     ..
1439 *     .. Local Scalars ..
1440       INTEGER            INTA, INTB, ZCODE
1441 *     ..
1442 *     .. Executable Statements ..
1443 *
1444 *     Test if the characters are equal
1445 *
1446       LSAME = CA.EQ.CB
1447       IF( LSAME )
1448      $   RETURN
1449 *
1450 *     Now test for equivalence if both characters are alphabetic.
1451 *
1452       ZCODE = ICHAR( 'Z' )
1453 *
1454 *     Use 'Z' rather than 'A' so that ASCII can be detected on Prime
1455 *     machines, on which ICHAR returns a value with bit 8 set.
1456 *     ICHAR('A') on Prime machines returns 193 which is the same as
1457 *     ICHAR('A') on an EBCDIC machine.
1458 *
1459       INTA = ICHAR( CA )
1460       INTB = ICHAR( CB )
1461 *
1462       IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN
1463 *
1464 *        ASCII is assumed - ZCODE is the ASCII code of either lower or
1465 *        upper case 'Z'.
1466 *
1467          IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32
1468          IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32
1469 *
1470       ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN
1471 *
1472 *        EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
1473 *        upper case 'Z'.
1474 *
1475          IF( INTA.GE.129 .AND. INTA.LE.137 .OR.
1476      $       INTA.GE.145 .AND. INTA.LE.153 .OR.
1477      $       INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64
1478          IF( INTB.GE.129 .AND. INTB.LE.137 .OR.
1479      $       INTB.GE.145 .AND. INTB.LE.153 .OR.
1480      $       INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64
1481 *
1482       ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN
1483 *
1484 *        ASCII is assumed, on Prime machines - ZCODE is the ASCII code
1485 *        plus 128 of either lower or upper case 'Z'.
1486 *
1487          IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32
1488          IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32
1489       END IF
1490       LSAME = INTA.EQ.INTB
1491 *
1492 *     RETURN
1493 *
1494 *     End of LSAME
1495 *
1496       END
1497       SUBROUTINE XERBLA( SRNAME, INFO )
1498 *
1499 *  -- LAPACK auxiliary routine (version 2.0) --
1500 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1501 *     Courant Institute, Argonne National Lab, and Rice University
1502 *     September 30, 1994
1503 *
1504 *     .. Scalar Arguments ..
1505       CHARACTER*6        SRNAME
1506       INTEGER            INFO
1507 *     ..
1508 *
1509 *  Purpose
1510 *  =======
1511 *
1512 *  XERBLA  is an error handler for the LAPACK routines.
1513 *  It is called by an LAPACK routine if an input parameter has an
1514 *  invalid value.  A message is printed and execution stops.
1515 *
1516 *  Installers may consider modifying the STOP statement in order to
1517 *  call system-specific exception-handling facilities.
1518 *
1519 *  Arguments
1520 *  =========
1521 *
1522 *  SRNAME  (input) CHARACTER*6
1523 *          The name of the routine which called XERBLA.
1524 *
1525 *  INFO    (input) INTEGER
1526 *          The position of the invalid parameter in the parameter list
1527 *          of the calling routine.
1528 *
1529 * =====================================================================
1530 *
1531 *     .. Executable Statements ..
1532 *
1533       WRITE( *, FMT = 9999 )SRNAME, INFO
1534 *
1535       STOP
1536 *
1537  9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ',
1538      $      'an illegal value' )
1539 *
1540 *     End of XERBLA
1541 *
1542       END
1543       SUBROUTINE DPOTF2( UPLO, N, A, LDA, INFO )
1544 *
1545 *  -- LAPACK routine (version 2.0) --
1546 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1547 *     Courant Institute, Argonne National Lab, and Rice University
1548 *     February 29, 1992
1549 *
1550 *     .. Scalar Arguments ..
1551       CHARACTER          UPLO
1552       INTEGER            INFO, LDA, N
1553 *     ..
1554 *     .. Array Arguments ..
1555       DOUBLE PRECISION   A( LDA, * )
1556 *     ..
1557 *
1558 *  Purpose
1559 *  =======
1560 *
1561 *  DPOTF2 computes the Cholesky factorization of a real symmetric
1562 *  positive definite matrix A.
1563 *
1564 *  The factorization has the form
1565 *     A = U' * U ,  if UPLO = 'U', or
1566 *     A = L  * L',  if UPLO = 'L',
1567 *  where U is an upper triangular matrix and L is lower triangular.
1568 *
1569 *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
1570 *
1571 *  Arguments
1572 *  =========
1573 *
1574 *  UPLO    (input) CHARACTER*1
1575 *          Specifies whether the upper or lower triangular part of the
1576 *          symmetric matrix A is stored.
1577 *          = 'U':  Upper triangular
1578 *          = 'L':  Lower triangular
1579 *
1580 *  N       (input) INTEGER
1581 *          The order of the matrix A.  N >= 0.
1582 *
1583 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
1584 *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
1585 *          n by n upper triangular part of A contains the upper
1586 *          triangular part of the matrix A, and the strictly lower
1587 *          triangular part of A is not referenced.  If UPLO = 'L', the
1588 *          leading n by n lower triangular part of A contains the lower
1589 *          triangular part of the matrix A, and the strictly upper
1590 *          triangular part of A is not referenced.
1591 *
1592 *          On exit, if INFO = 0, the factor U or L from the Cholesky
1593 *          factorization A = U'*U  or A = L*L'.
1594 *
1595 *  LDA     (input) INTEGER
1596 *          The leading dimension of the array A.  LDA >= max(1,N).
1597 *
1598 *  INFO    (output) INTEGER
1599 *          = 0: successful exit
1600 *          < 0: if INFO = -k, the k-th argument had an illegal value
1601 *          > 0: if INFO = k, the leading minor of order k is not
1602 *               positive definite, and the factorization could not be
1603 *               completed.
1604 *
1605 *  =====================================================================
1606 *
1607 *     .. Parameters ..
1608       DOUBLE PRECISION   ONE, ZERO
1609       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
1610 *     ..
1611 *     .. Local Scalars ..
1612       LOGICAL            UPPER
1613       INTEGER            J
1614       DOUBLE PRECISION   AJJ
1615 *     ..
1616 *     .. External Functions ..
1617       LOGICAL            LSAME
1618       DOUBLE PRECISION   DDOT
1619       EXTERNAL           LSAME, DDOT
1620 *     ..
1621 *     .. External Subroutines ..
1622       EXTERNAL           DGEMV, DSCAL, XERBLA
1623 *     ..
1624 *     .. Intrinsic Functions ..
1625       INTRINSIC          MAX, SQRT
1626 *     ..
1627 *     .. Executable Statements ..
1628 *
1629 *     Test the input parameters.
1630 *
1631       INFO = 0
1632       UPPER = LSAME( UPLO, 'U' )
1633       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
1634          INFO = -1
1635       ELSE IF( N.LT.0 ) THEN
1636          INFO = -2
1637       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
1638          INFO = -4
1639       END IF
1640       IF( INFO.NE.0 ) THEN
1641          CALL XERBLA( 'DPOTF2', -INFO )
1642          RETURN
1643       END IF
1644 *
1645 *     Quick return if possible
1646 *
1647       IF( N.EQ.0 )
1648      $   RETURN
1649 *
1650       IF( UPPER ) THEN
1651 *
1652 *        Compute the Cholesky factorization A = U'*U.
1653 *
1654          DO 10 J = 1, N
1655 *
1656 *           Compute U(J,J) and test for non-positive-definiteness.
1657 *
1658             AJJ = A( J, J ) - DDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 )
1659             IF( AJJ.LE.ZERO ) THEN
1660                A( J, J ) = AJJ
1661                GO TO 30
1662             END IF
1663             AJJ = SQRT( AJJ )
1664             A( J, J ) = AJJ
1665 *
1666 *           Compute elements J+1:N of row J.
1667 *
1668             IF( J.LT.N ) THEN
1669                CALL DGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ),
1670      $                     LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
1671                CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
1672             END IF
1673    10    CONTINUE
1674       ELSE
1675 *
1676 *        Compute the Cholesky factorization A = L*L'.
1677 *
1678          DO 20 J = 1, N
1679 *
1680 *           Compute L(J,J) and test for non-positive-definiteness.
1681 *
1682             AJJ = A( J, J ) - DDOT( J-1, A( J, 1 ), LDA, A( J, 1 ),
1683      $            LDA )
1684             IF( AJJ.LE.ZERO ) THEN
1685                A( J, J ) = AJJ
1686                GO TO 30
1687             END IF
1688             AJJ = SQRT( AJJ )
1689             A( J, J ) = AJJ
1690 *
1691 *           Compute elements J+1:N of column J.
1692 *
1693             IF( J.LT.N ) THEN
1694                CALL DGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ),
1695      $                     LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
1696                CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
1697             END IF
1698    20    CONTINUE
1699       END IF
1700       GO TO 40
1701 *
1702    30 CONTINUE
1703       INFO = J
1704 *
1705    40 CONTINUE
1706       RETURN
1707 *
1708 *     End of DPOTF2
1709 *
1710       END
1711       subroutine  dscal(n,da,dx,incx)
1712 c
1713 c     scales a vector by a constant.
1714 c     uses unrolled loops for increment equal to one.
1715 c     jack dongarra, linpack, 3/11/78.
1716 c     modified 3/93 to return if incx .le. 0.
1717 c     modified 12/3/93, array(1) declarations changed to array(*)
1718 c
1719       double precision da,dx(*)
1720       integer i,incx,m,mp1,n,nincx
1721 c
1722       if( n.le.0 .or. incx.le.0 )return
1723       if(incx.eq.1)go to 20
1724 c
1725 c        code for increment not equal to 1
1726 c
1727       nincx = n*incx
1728       do 10 i = 1,nincx,incx
1729         dx(i) = da*dx(i)
1730    10 continue
1731       return
1732 c
1733 c        code for increment equal to 1
1734 c
1735 c
1736 c        clean-up loop
1737 c
1738    20 m = mod(n,5)
1739       if( m .eq. 0 ) go to 40
1740       do 30 i = 1,m
1741         dx(i) = da*dx(i)
1742    30 continue
1743       if( n .lt. 5 ) return
1744    40 mp1 = m + 1
1745       do 50 i = mp1,n,5
1746         dx(i) = da*dx(i)
1747         dx(i + 1) = da*dx(i + 1)
1748         dx(i + 2) = da*dx(i + 2)
1749         dx(i + 3) = da*dx(i + 3)
1750         dx(i + 4) = da*dx(i + 4)
1751    50 continue
1752       return
1753       end
1754       double precision function ddot(n,dx,incx,dy,incy)
1755 c
1756 c     forms the dot product of two vectors.
1757 c     uses unrolled loops for increments equal to one.
1758 c     jack dongarra, linpack, 3/11/78.
1759 c     modified 12/3/93, array(1) declarations changed to array(*)
1760 c
1761       double precision dx(*),dy(*),dtemp
1762       integer i,incx,incy,ix,iy,m,mp1,n
1763 c
1764       ddot = 0.0d0
1765       dtemp = 0.0d0
1766       if(n.le.0)return
1767       if(incx.eq.1.and.incy.eq.1)go to 20
1768 c
1769 c        code for unequal increments or equal increments
1770 c          not equal to 1
1771 c
1772       ix = 1
1773       iy = 1
1774       if(incx.lt.0)ix = (-n+1)*incx + 1
1775       if(incy.lt.0)iy = (-n+1)*incy + 1
1776       do 10 i = 1,n
1777         dtemp = dtemp + dx(ix)*dy(iy)
1778         ix = ix + incx
1779         iy = iy + incy
1780    10 continue
1781       ddot = dtemp
1782       return
1783 c
1784 c        code for both increments equal to 1
1785 c
1786 c
1787 c        clean-up loop
1788 c
1789    20 m = mod(n,5)
1790       if( m .eq. 0 ) go to 40
1791       do 30 i = 1,m
1792         dtemp = dtemp + dx(i)*dy(i)
1793    30 continue
1794       if( n .lt. 5 ) go to 60
1795    40 mp1 = m + 1
1796       do 50 i = mp1,n,5
1797         dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
1798      *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
1799    50 continue
1800    60 ddot = dtemp
1801       return
1802       end
1803       SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
1804      $                   B, LDB )
1805 *     .. Scalar Arguments ..
1806       CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
1807       INTEGER            M, N, LDA, LDB
1808       DOUBLE PRECISION   ALPHA
1809 *     .. Array Arguments ..
1810       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
1811 *     ..
1812 *
1813 *  Purpose
1814 *  =======
1815 *
1816 *  DTRSM  solves one of the matrix equations
1817 *
1818 *     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
1819 *
1820 *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
1821 *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
1822 *
1823 *     op( A ) = A   or   op( A ) = A'.
1824 *
1825 *  The matrix X is overwritten on B.
1826 *
1827 *  Parameters
1828 *  ==========
1829 *
1830 *  SIDE   - CHARACTER*1.
1831 *           On entry, SIDE specifies whether op( A ) appears on the left
1832 *           or right of X as follows:
1833 *
1834 *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
1835 *
1836 *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
1837 *
1838 *           Unchanged on exit.
1839 *
1840 *  UPLO   - CHARACTER*1.
1841 *           On entry, UPLO specifies whether the matrix A is an upper or
1842 *           lower triangular matrix as follows:
1843 *
1844 *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
1845 *
1846 *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
1847 *
1848 *           Unchanged on exit.
1849 *
1850 *  TRANSA - CHARACTER*1.
1851 *           On entry, TRANSA specifies the form of op( A ) to be used in
1852 *           the matrix multiplication as follows:
1853 *
1854 *              TRANSA = 'N' or 'n'   op( A ) = A.
1855 *
1856 *              TRANSA = 'T' or 't'   op( A ) = A'.
1857 *
1858 *              TRANSA = 'C' or 'c'   op( A ) = A'.
1859 *
1860 *           Unchanged on exit.
1861 *
1862 *  DIAG   - CHARACTER*1.
1863 *           On entry, DIAG specifies whether or not A is unit triangular
1864 *           as follows:
1865 *
1866 *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
1867 *
1868 *              DIAG = 'N' or 'n'   A is not assumed to be unit
1869 *                                  triangular.
1870 *
1871 *           Unchanged on exit.
1872 *
1873 *  M      - INTEGER.
1874 *           On entry, M specifies the number of rows of B. M must be at
1875 *           least zero.
1876 *           Unchanged on exit.
1877 *
1878 *  N      - INTEGER.
1879 *           On entry, N specifies the number of columns of B.  N must be
1880 *           at least zero.
1881 *           Unchanged on exit.
1882 *
1883 *  ALPHA  - DOUBLE PRECISION.
1884 *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
1885 *           zero then  A is not referenced and  B need not be set before
1886 *           entry.
1887 *           Unchanged on exit.
1888 *
1889 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
1890 *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
1891 *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
1892 *           upper triangular part of the array  A must contain the upper
1893 *           triangular matrix  and the strictly lower triangular part of
1894 *           A is not referenced.
1895 *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
1896 *           lower triangular part of the array  A must contain the lower
1897 *           triangular matrix  and the strictly upper triangular part of
1898 *           A is not referenced.
1899 *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
1900 *           A  are not referenced either,  but are assumed to be  unity.
1901 *           Unchanged on exit.
1902 *
1903 *  LDA    - INTEGER.
1904 *           On entry, LDA specifies the first dimension of A as declared
1905 *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
1906 *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
1907 *           then LDA must be at least max( 1, n ).
1908 *           Unchanged on exit.
1909 *
1910 *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
1911 *           Before entry,  the leading  m by n part of the array  B must
1912 *           contain  the  right-hand  side  matrix  B,  and  on exit  is
1913 *           overwritten by the solution matrix  X.
1914 *
1915 *  LDB    - INTEGER.
1916 *           On entry, LDB specifies the first dimension of B as declared
1917 *           in  the  calling  (sub)  program.   LDB  must  be  at  least
1918 *           max( 1, m ).
1919 *           Unchanged on exit.
1920 *
1921 *
1922 *  Level 3 Blas routine.
1923 *
1924 *
1925 *  -- Written on 8-February-1989.
1926 *     Jack Dongarra, Argonne National Laboratory.
1927 *     Iain Duff, AERE Harwell.
1928 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
1929 *     Sven Hammarling, Numerical Algorithms Group Ltd.
1930 *
1931 *
1932 *     .. External Functions ..
1933       LOGICAL            LSAME
1934       EXTERNAL           LSAME
1935 *     .. External Subroutines ..
1936       EXTERNAL           XERBLA
1937 *     .. Intrinsic Functions ..
1938       INTRINSIC          MAX
1939 *     .. Local Scalars ..
1940       LOGICAL            LSIDE, NOUNIT, UPPER
1941       INTEGER            I, INFO, J, K, NROWA
1942       DOUBLE PRECISION   TEMP
1943 *     .. Parameters ..
1944       DOUBLE PRECISION   ONE         , ZERO
1945       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
1946 *     ..
1947 *     .. Executable Statements ..
1948 *
1949 *     Test the input parameters.
1950 *
1951       LSIDE  = LSAME( SIDE  , 'L' )
1952       IF( LSIDE )THEN
1953          NROWA = M
1954       ELSE
1955          NROWA = N
1956       END IF
1957       NOUNIT = LSAME( DIAG  , 'N' )
1958       UPPER  = LSAME( UPLO  , 'U' )
1959 *
1960       INFO   = 0
1961       IF(      ( .NOT.LSIDE                ).AND.
1962      $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
1963          INFO = 1
1964       ELSE IF( ( .NOT.UPPER                ).AND.
1965      $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
1966          INFO = 2
1967       ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
1968      $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
1969      $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
1970          INFO = 3
1971       ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
1972      $         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
1973          INFO = 4
1974       ELSE IF( M  .LT.0               )THEN
1975          INFO = 5
1976       ELSE IF( N  .LT.0               )THEN
1977          INFO = 6
1978       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
1979          INFO = 9
1980       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
1981          INFO = 11
1982       END IF
1983       IF( INFO.NE.0 )THEN
1984          CALL XERBLA( 'DTRSM ', INFO )
1985          RETURN
1986       END IF
1987 *
1988 *     Quick return if possible.
1989 *
1990       IF( N.EQ.0 )
1991      $   RETURN
1992 *
1993 *     And when  alpha.eq.zero.
1994 *
1995       IF( ALPHA.EQ.ZERO )THEN
1996          DO 20, J = 1, N
1997             DO 10, I = 1, M
1998                B( I, J ) = ZERO
1999    10       CONTINUE
2000    20    CONTINUE
2001          RETURN
2002       END IF
2003 *
2004 *     Start the operations.
2005 *
2006       IF( LSIDE )THEN
2007          IF( LSAME( TRANSA, 'N' ) )THEN
2008 *
2009 *           Form  B := alpha*inv( A )*B.
2010 *
2011             IF( UPPER )THEN
2012                DO 60, J = 1, N
2013                   IF( ALPHA.NE.ONE )THEN
2014                      DO 30, I = 1, M
2015                         B( I, J ) = ALPHA*B( I, J )
2016    30                CONTINUE
2017                   END IF
2018                   DO 50, K = M, 1, -1
2019                      IF( B( K, J ).NE.ZERO )THEN
2020                         IF( NOUNIT )
2021      $                     B( K, J ) = B( K, J )/A( K, K )
2022                         DO 40, I = 1, K - 1
2023                            B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
2024    40                   CONTINUE
2025                      END IF
2026    50             CONTINUE
2027    60          CONTINUE
2028             ELSE
2029                DO 100, J = 1, N
2030                   IF( ALPHA.NE.ONE )THEN
2031                      DO 70, I = 1, M
2032                         B( I, J ) = ALPHA*B( I, J )
2033    70                CONTINUE
2034                   END IF
2035                   DO 90 K = 1, M
2036                      IF( B( K, J ).NE.ZERO )THEN
2037                         IF( NOUNIT )
2038      $                     B( K, J ) = B( K, J )/A( K, K )
2039                         DO 80, I = K + 1, M
2040                            B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
2041    80                   CONTINUE
2042                      END IF
2043    90             CONTINUE
2044   100          CONTINUE
2045             END IF
2046          ELSE
2047 *
2048 *           Form  B := alpha*inv( A' )*B.
2049 *
2050             IF( UPPER )THEN
2051                DO 130, J = 1, N
2052                   DO 120, I = 1, M
2053                      TEMP = ALPHA*B( I, J )
2054                      DO 110, K = 1, I - 1
2055                         TEMP = TEMP - A( K, I )*B( K, J )
2056   110                CONTINUE
2057                      IF( NOUNIT )
2058      $                  TEMP = TEMP/A( I, I )
2059                      B( I, J ) = TEMP
2060   120             CONTINUE
2061   130          CONTINUE
2062             ELSE
2063                DO 160, J = 1, N
2064                   DO 150, I = M, 1, -1
2065                      TEMP = ALPHA*B( I, J )
2066                      DO 140, K = I + 1, M
2067                         TEMP = TEMP - A( K, I )*B( K, J )
2068   140                CONTINUE
2069                      IF( NOUNIT )
2070      $                  TEMP = TEMP/A( I, I )
2071                      B( I, J ) = TEMP
2072   150             CONTINUE
2073   160          CONTINUE
2074             END IF
2075          END IF
2076       ELSE
2077          IF( LSAME( TRANSA, 'N' ) )THEN
2078 *
2079 *           Form  B := alpha*B*inv( A ).
2080 *
2081             IF( UPPER )THEN
2082                DO 210, J = 1, N
2083                   IF( ALPHA.NE.ONE )THEN
2084                      DO 170, I = 1, M
2085                         B( I, J ) = ALPHA*B( I, J )
2086   170                CONTINUE
2087                   END IF
2088                   DO 190, K = 1, J - 1
2089                      IF( A( K, J ).NE.ZERO )THEN
2090                         DO 180, I = 1, M
2091                            B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
2092   180                   CONTINUE
2093                      END IF
2094   190             CONTINUE
2095                   IF( NOUNIT )THEN
2096                      TEMP = ONE/A( J, J )
2097                      DO 200, I = 1, M
2098                         B( I, J ) = TEMP*B( I, J )
2099   200                CONTINUE
2100                   END IF
2101   210          CONTINUE
2102             ELSE
2103                DO 260, J = N, 1, -1
2104                   IF( ALPHA.NE.ONE )THEN
2105                      DO 220, I = 1, M
2106                         B( I, J ) = ALPHA*B( I, J )
2107   220                CONTINUE
2108                   END IF
2109                   DO 240, K = J + 1, N
2110                      IF( A( K, J ).NE.ZERO )THEN
2111                         DO 230, I = 1, M
2112                            B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
2113   230                   CONTINUE
2114                      END IF
2115   240             CONTINUE
2116                   IF( NOUNIT )THEN
2117                      TEMP = ONE/A( J, J )
2118                      DO 250, I = 1, M
2119                        B( I, J ) = TEMP*B( I, J )
2120   250                CONTINUE
2121                   END IF
2122   260          CONTINUE
2123             END IF
2124          ELSE
2125 *
2126 *           Form  B := alpha*B*inv( A' ).
2127 *
2128             IF( UPPER )THEN
2129                DO 310, K = N, 1, -1
2130                   IF( NOUNIT )THEN
2131                      TEMP = ONE/A( K, K )
2132                      DO 270, I = 1, M
2133                         B( I, K ) = TEMP*B( I, K )
2134   270                CONTINUE
2135                   END IF
2136                   DO 290, J = 1, K - 1
2137                      IF( A( J, K ).NE.ZERO )THEN
2138                         TEMP = A( J, K )
2139                         DO 280, I = 1, M
2140                            B( I, J ) = B( I, J ) - TEMP*B( I, K )
2141   280                   CONTINUE
2142                      END IF
2143   290             CONTINUE
2144                   IF( ALPHA.NE.ONE )THEN
2145                      DO 300, I = 1, M
2146                         B( I, K ) = ALPHA*B( I, K )
2147   300                CONTINUE
2148                   END IF
2149   310          CONTINUE
2150             ELSE
2151                DO 360, K = 1, N
2152                   IF( NOUNIT )THEN
2153                      TEMP = ONE/A( K, K )
2154                      DO 320, I = 1, M
2155                         B( I, K ) = TEMP*B( I, K )
2156   320                CONTINUE
2157                   END IF
2158                   DO 340, J = K + 1, N
2159                      IF( A( J, K ).NE.ZERO )THEN
2160                         TEMP = A( J, K )
2161                         DO 330, I = 1, M
2162                            B( I, J ) = B( I, J ) - TEMP*B( I, K )
2163   330                   CONTINUE
2164                      END IF
2165   340             CONTINUE
2166                   IF( ALPHA.NE.ONE )THEN
2167                      DO 350, I = 1, M
2168                         B( I, K ) = ALPHA*B( I, K )
2169   350                CONTINUE
2170                   END IF
2171   360          CONTINUE
2172             END IF
2173          END IF
2174       END IF
2175 *
2176       RETURN
2177 *
2178 *     End of DTRSM .
2179 *
2180       END
2181       SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
2182 *     .. Scalar Arguments ..
2183       INTEGER            INCX, K, LDA, N
2184       CHARACTER*1        DIAG, TRANS, UPLO
2185 *     .. Array Arguments ..
2186       DOUBLE PRECISION   A( LDA, * ), X( * )
2187 *     ..
2188 *
2189 *  Purpose
2190 *  =======
2191 *
2192 *  DTBSV  solves one of the systems of equations
2193 *
2194 *     A*x = b,   or   A'*x = b,
2195 *
2196 *  where b and x are n element vectors and A is an n by n unit, or
2197 *  non-unit, upper or lower triangular band matrix, with ( k + 1 )
2198 *  diagonals.
2199 *
2200 *  No test for singularity or near-singularity is included in this
2201 *  routine. Such tests must be performed before calling this routine.
2202 *
2203 *  Parameters
2204 *  ==========
2205 *
2206 *  UPLO   - CHARACTER*1.
2207 *           On entry, UPLO specifies whether the matrix is an upper or
2208 *           lower triangular matrix as follows:
2209 *
2210 *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
2211 *
2212 *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
2213 *
2214 *           Unchanged on exit.
2215 *
2216 *  TRANS  - CHARACTER*1.
2217 *           On entry, TRANS specifies the equations to be solved as
2218 *           follows:
2219 *
2220 *              TRANS = 'N' or 'n'   A*x = b.
2221 *
2222 *              TRANS = 'T' or 't'   A'*x = b.
2223 *
2224 *              TRANS = 'C' or 'c'   A'*x = b.
2225 *
2226 *           Unchanged on exit.
2227 *
2228 *  DIAG   - CHARACTER*1.
2229 *           On entry, DIAG specifies whether or not A is unit
2230 *           triangular as follows:
2231 *
2232 *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
2233 *
2234 *              DIAG = 'N' or 'n'   A is not assumed to be unit
2235 *                                  triangular.
2236 *
2237 *           Unchanged on exit.
2238 *
2239 *  N      - INTEGER.
2240 *           On entry, N specifies the order of the matrix A.
2241 *           N must be at least zero.
2242 *           Unchanged on exit.
2243 *
2244 *  K      - INTEGER.
2245 *           On entry with UPLO = 'U' or 'u', K specifies the number of
2246 *           super-diagonals of the matrix A.
2247 *           On entry with UPLO = 'L' or 'l', K specifies the number of
2248 *           sub-diagonals of the matrix A.
2249 *           K must satisfy  0 .le. K.
2250 *           Unchanged on exit.
2251 *
2252 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2253 *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
2254 *           by n part of the array A must contain the upper triangular
2255 *           band part of the matrix of coefficients, supplied column by
2256 *           column, with the leading diagonal of the matrix in row
2257 *           ( k + 1 ) of the array, the first super-diagonal starting at
2258 *           position 2 in row k, and so on. The top left k by k triangle
2259 *           of the array A is not referenced.
2260 *           The following program segment will transfer an upper
2261 *           triangular band matrix from conventional full matrix storage
2262 *           to band storage:
2263 *
2264 *                 DO 20, J = 1, N
2265 *                    M = K + 1 - J
2266 *                    DO 10, I = MAX( 1, J - K ), J
2267 *                       A( M + I, J ) = matrix( I, J )
2268 *              10    CONTINUE
2269 *              20 CONTINUE
2270 *
2271 *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
2272 *           by n part of the array A must contain the lower triangular
2273 *           band part of the matrix of coefficients, supplied column by
2274 *           column, with the leading diagonal of the matrix in row 1 of
2275 *           the array, the first sub-diagonal starting at position 1 in
2276 *           row 2, and so on. The bottom right k by k triangle of the
2277 *           array A is not referenced.
2278 *           The following program segment will transfer a lower
2279 *           triangular band matrix from conventional full matrix storage
2280 *           to band storage:
2281 *
2282 *                 DO 20, J = 1, N
2283 *                    M = 1 - J
2284 *                    DO 10, I = J, MIN( N, J + K )
2285 *                       A( M + I, J ) = matrix( I, J )
2286 *              10    CONTINUE
2287 *              20 CONTINUE
2288 *
2289 *           Note that when DIAG = 'U' or 'u' the elements of the array A
2290 *           corresponding to the diagonal elements of the matrix are not
2291 *           referenced, but are assumed to be unity.
2292 *           Unchanged on exit.
2293 *
2294 *  LDA    - INTEGER.
2295 *           On entry, LDA specifies the first dimension of A as declared
2296 *           in the calling (sub) program. LDA must be at least
2297 *           ( k + 1 ).
2298 *           Unchanged on exit.
2299 *
2300 *  X      - DOUBLE PRECISION array of dimension at least
2301 *           ( 1 + ( n - 1 )*abs( INCX ) ).
2302 *           Before entry, the incremented array X must contain the n
2303 *           element right-hand side vector b. On exit, X is overwritten
2304 *           with the solution vector x.
2305 *
2306 *  INCX   - INTEGER.
2307 *           On entry, INCX specifies the increment for the elements of
2308 *           X. INCX must not be zero.
2309 *           Unchanged on exit.
2310 *
2311 *
2312 *  Level 2 Blas routine.
2313 *
2314 *  -- Written on 22-October-1986.
2315 *     Jack Dongarra, Argonne National Lab.
2316 *     Jeremy Du Croz, Nag Central Office.
2317 *     Sven Hammarling, Nag Central Office.
2318 *     Richard Hanson, Sandia National Labs.
2319 *
2320 *
2321 *     .. Parameters ..
2322       DOUBLE PRECISION   ZERO
2323       PARAMETER        ( ZERO = 0.0D+0 )
2324 *     .. Local Scalars ..
2325       DOUBLE PRECISION   TEMP
2326       INTEGER            I, INFO, IX, J, JX, KPLUS1, KX, L
2327       LOGICAL            NOUNIT
2328 *     .. External Functions ..
2329       LOGICAL            LSAME
2330       EXTERNAL           LSAME
2331 *     .. External Subroutines ..
2332       EXTERNAL           XERBLA
2333 *     .. Intrinsic Functions ..
2334       INTRINSIC          MAX, MIN
2335 *     ..
2336 *     .. Executable Statements ..
2337 *
2338 *     Test the input parameters.
2339 *
2340       INFO = 0
2341       IF     ( .NOT.LSAME( UPLO , 'U' ).AND.
2342      $         .NOT.LSAME( UPLO , 'L' )      )THEN
2343          INFO = 1
2344       ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
2345      $         .NOT.LSAME( TRANS, 'T' ).AND.
2346      $         .NOT.LSAME( TRANS, 'C' )      )THEN
2347          INFO = 2
2348       ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
2349      $         .NOT.LSAME( DIAG , 'N' )      )THEN
2350          INFO = 3
2351       ELSE IF( N.LT.0 )THEN
2352          INFO = 4
2353       ELSE IF( K.LT.0 )THEN
2354          INFO = 5
2355       ELSE IF( LDA.LT.( K + 1 ) )THEN
2356          INFO = 7
2357       ELSE IF( INCX.EQ.0 )THEN
2358          INFO = 9
2359       END IF
2360       IF( INFO.NE.0 )THEN
2361          CALL XERBLA( 'DTBSV ', INFO )
2362          RETURN
2363       END IF
2364 *
2365 *     Quick return if possible.
2366 *
2367       IF( N.EQ.0 )
2368      $   RETURN
2369 *
2370       NOUNIT = LSAME( DIAG, 'N' )
2371 *
2372 *     Set up the start point in X if the increment is not unity. This
2373 *     will be  ( N - 1 )*INCX  too small for descending loops.
2374 *
2375       IF( INCX.LE.0 )THEN
2376          KX = 1 - ( N - 1 )*INCX
2377       ELSE IF( INCX.NE.1 )THEN
2378          KX = 1
2379       END IF
2380 *
2381 *     Start the operations. In this version the elements of A are
2382 *     accessed by sequentially with one pass through A.
2383 *
2384       IF( LSAME( TRANS, 'N' ) )THEN
2385 *
2386 *        Form  x := inv( A )*x.
2387 *
2388          IF( LSAME( UPLO, 'U' ) )THEN
2389             KPLUS1 = K + 1
2390             IF( INCX.EQ.1 )THEN
2391                DO 20, J = N, 1, -1
2392                   IF( X( J ).NE.ZERO )THEN
2393                      L = KPLUS1 - J
2394                      IF( NOUNIT )
2395      $                  X( J ) = X( J )/A( KPLUS1, J )
2396                      TEMP = X( J )
2397                      DO 10, I = J - 1, MAX( 1, J - K ), -1
2398                         X( I ) = X( I ) - TEMP*A( L + I, J )
2399    10                CONTINUE
2400                   END IF
2401    20          CONTINUE
2402             ELSE
2403                KX = KX + ( N - 1 )*INCX
2404                JX = KX
2405                DO 40, J = N, 1, -1
2406                   KX = KX - INCX
2407                   IF( X( JX ).NE.ZERO )THEN
2408                      IX = KX
2409                      L  = KPLUS1 - J
2410                      IF( NOUNIT )
2411      $                  X( JX ) = X( JX )/A( KPLUS1, J )
2412                      TEMP = X( JX )
2413                      DO 30, I = J - 1, MAX( 1, J - K ), -1
2414                         X( IX ) = X( IX ) - TEMP*A( L + I, J )
2415                         IX      = IX      - INCX
2416    30                CONTINUE
2417                   END IF
2418                   JX = JX - INCX
2419    40          CONTINUE
2420             END IF
2421          ELSE
2422             IF( INCX.EQ.1 )THEN
2423                DO 60, J = 1, N
2424                   IF( X( J ).NE.ZERO )THEN
2425                      L = 1 - J
2426                      IF( NOUNIT )
2427      $                  X( J ) = X( J )/A( 1, J )
2428                      TEMP = X( J )
2429                      DO 50, I = J + 1, MIN( N, J + K )
2430                         X( I ) = X( I ) - TEMP*A( L + I, J )
2431    50                CONTINUE
2432                   END IF
2433    60          CONTINUE
2434             ELSE
2435                JX = KX
2436                DO 80, J = 1, N
2437                   KX = KX + INCX
2438                   IF( X( JX ).NE.ZERO )THEN
2439                      IX = KX
2440                      L  = 1  - J
2441                      IF( NOUNIT )
2442      $                  X( JX ) = X( JX )/A( 1, J )
2443                      TEMP = X( JX )
2444                      DO 70, I = J + 1, MIN( N, J + K )
2445                         X( IX ) = X( IX ) - TEMP*A( L + I, J )
2446                         IX      = IX      + INCX
2447    70                CONTINUE
2448                   END IF
2449                   JX = JX + INCX
2450    80          CONTINUE
2451             END IF
2452          END IF
2453       ELSE
2454 *
2455 *        Form  x := inv( A')*x.
2456 *
2457          IF( LSAME( UPLO, 'U' ) )THEN
2458             KPLUS1 = K + 1
2459             IF( INCX.EQ.1 )THEN
2460                DO 100, J = 1, N
2461                   TEMP = X( J )
2462                   L    = KPLUS1 - J
2463                   DO 90, I = MAX( 1, J - K ), J - 1
2464                      TEMP = TEMP - A( L + I, J )*X( I )
2465    90             CONTINUE
2466                   IF( NOUNIT )
2467      $               TEMP = TEMP/A( KPLUS1, J )
2468                   X( J ) = TEMP
2469   100          CONTINUE
2470             ELSE
2471                JX = KX
2472                DO 120, J = 1, N
2473                   TEMP = X( JX )
2474                   IX   = KX
2475                   L    = KPLUS1  - J
2476                   DO 110, I = MAX( 1, J - K ), J - 1
2477                      TEMP = TEMP - A( L + I, J )*X( IX )
2478                      IX   = IX   + INCX
2479   110             CONTINUE
2480                   IF( NOUNIT )
2481      $               TEMP = TEMP/A( KPLUS1, J )
2482                   X( JX ) = TEMP
2483                   JX      = JX   + INCX
2484                   IF( J.GT.K )
2485      $               KX = KX + INCX
2486   120          CONTINUE
2487             END IF
2488          ELSE
2489             IF( INCX.EQ.1 )THEN
2490                DO 140, J = N, 1, -1
2491                   TEMP = X( J )
2492                   L    = 1      - J
2493                   DO 130, I = MIN( N, J + K ), J + 1, -1
2494                      TEMP = TEMP - A( L + I, J )*X( I )
2495   130             CONTINUE
2496                   IF( NOUNIT )
2497      $               TEMP = TEMP/A( 1, J )
2498                   X( J ) = TEMP
2499   140          CONTINUE
2500             ELSE
2501                KX = KX + ( N - 1 )*INCX
2502                JX = KX
2503                DO 160, J = N, 1, -1
2504                   TEMP = X( JX )
2505                   IX   = KX
2506                   L    = 1       - J
2507                   DO 150, I = MIN( N, J + K ), J + 1, -1
2508                      TEMP = TEMP - A( L + I, J )*X( IX )
2509                      IX   = IX   - INCX
2510   150             CONTINUE
2511                   IF( NOUNIT )
2512      $               TEMP = TEMP/A( 1, J )
2513                   X( JX ) = TEMP
2514                   JX      = JX   - INCX
2515                   IF( ( N - J ).GE.K )
2516      $               KX = KX - INCX
2517   160          CONTINUE
2518             END IF
2519          END IF
2520       END IF
2521 *
2522       RETURN
2523 *
2524 *     End of DTBSV .
2525 *
2526       END
2527       SUBROUTINE DSYR  ( UPLO, N, ALPHA, X, INCX, A, LDA )
2528 *     .. Scalar Arguments ..
2529       DOUBLE PRECISION   ALPHA
2530       INTEGER            INCX, LDA, N
2531       CHARACTER*1        UPLO
2532 *     .. Array Arguments ..
2533       DOUBLE PRECISION   A( LDA, * ), X( * )
2534 *     ..
2535 *
2536 *  Purpose
2537 *  =======
2538 *
2539 *  DSYR   performs the symmetric rank 1 operation
2540 *
2541 *     A := alpha*x*x' + A,
2542 *
2543 *  where alpha is a real scalar, x is an n element vector and A is an
2544 *  n by n symmetric matrix.
2545 *
2546 *  Parameters
2547 *  ==========
2548 *
2549 *  UPLO   - CHARACTER*1.
2550 *           On entry, UPLO specifies whether the upper or lower
2551 *           triangular part of the array A is to be referenced as
2552 *           follows:
2553 *
2554 *              UPLO = 'U' or 'u'   Only the upper triangular part of A
2555 *                                  is to be referenced.
2556 *
2557 *              UPLO = 'L' or 'l'   Only the lower triangular part of A
2558 *                                  is to be referenced.
2559 *
2560 *           Unchanged on exit.
2561 *
2562 *  N      - INTEGER.
2563 *           On entry, N specifies the order of the matrix A.
2564 *           N must be at least zero.
2565 *           Unchanged on exit.
2566 *
2567 *  ALPHA  - DOUBLE PRECISION.
2568 *           On entry, ALPHA specifies the scalar alpha.
2569 *           Unchanged on exit.
2570 *
2571 *  X      - DOUBLE PRECISION array of dimension at least
2572 *           ( 1 + ( n - 1 )*abs( INCX ) ).
2573 *           Before entry, the incremented array X must contain the n
2574 *           element vector x.
2575 *           Unchanged on exit.
2576 *
2577 *  INCX   - INTEGER.
2578 *           On entry, INCX specifies the increment for the elements of
2579 *           X. INCX must not be zero.
2580 *           Unchanged on exit.
2581 *
2582 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2583 *           Before entry with  UPLO = 'U' or 'u', the leading n by n
2584 *           upper triangular part of the array A must contain the upper
2585 *           triangular part of the symmetric matrix and the strictly
2586 *           lower triangular part of A is not referenced. On exit, the
2587 *           upper triangular part of the array A is overwritten by the
2588 *           upper triangular part of the updated matrix.
2589 *           Before entry with UPLO = 'L' or 'l', the leading n by n
2590 *           lower triangular part of the array A must contain the lower
2591 *           triangular part of the symmetric matrix and the strictly
2592 *           upper triangular part of A is not referenced. On exit, the
2593 *           lower triangular part of the array A is overwritten by the
2594 *           lower triangular part of the updated matrix.
2595 *
2596 *  LDA    - INTEGER.
2597 *           On entry, LDA specifies the first dimension of A as declared
2598 *           in the calling (sub) program. LDA must be at least
2599 *           max( 1, n ).
2600 *           Unchanged on exit.
2601 *
2602 *
2603 *  Level 2 Blas routine.
2604 *
2605 *  -- Written on 22-October-1986.
2606 *     Jack Dongarra, Argonne National Lab.
2607 *     Jeremy Du Croz, Nag Central Office.
2608 *     Sven Hammarling, Nag Central Office.
2609 *     Richard Hanson, Sandia National Labs.
2610 *
2611 *
2612 *     .. Parameters ..
2613       DOUBLE PRECISION   ZERO
2614       PARAMETER        ( ZERO = 0.0D+0 )
2615 *     .. Local Scalars ..
2616       DOUBLE PRECISION   TEMP
2617       INTEGER            I, INFO, IX, J, JX, KX
2618 *     .. External Functions ..
2619       LOGICAL            LSAME
2620       EXTERNAL           LSAME
2621 *     .. External Subroutines ..
2622       EXTERNAL           XERBLA
2623 *     .. Intrinsic Functions ..
2624       INTRINSIC          MAX
2625 *     ..
2626 *     .. Executable Statements ..
2627 *
2628 *     Test the input parameters.
2629 *
2630       INFO = 0
2631       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
2632      $         .NOT.LSAME( UPLO, 'L' )      )THEN
2633          INFO = 1
2634       ELSE IF( N.LT.0 )THEN
2635          INFO = 2
2636       ELSE IF( INCX.EQ.0 )THEN
2637          INFO = 5
2638       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
2639          INFO = 7
2640       END IF
2641       IF( INFO.NE.0 )THEN
2642          CALL XERBLA( 'DSYR  ', INFO )
2643          RETURN
2644       END IF
2645 *
2646 *     Quick return if possible.
2647 *
2648       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
2649      $   RETURN
2650 *
2651 *     Set the start point in X if the increment is not unity.
2652 *
2653       IF( INCX.LE.0 )THEN
2654          KX = 1 - ( N - 1 )*INCX
2655       ELSE IF( INCX.NE.1 )THEN
2656          KX = 1
2657       END IF
2658 *
2659 *     Start the operations. In this version the elements of A are
2660 *     accessed sequentially with one pass through the triangular part
2661 *     of A.
2662 *
2663       IF( LSAME( UPLO, 'U' ) )THEN
2664 *
2665 *        Form  A  when A is stored in upper triangle.
2666 *
2667          IF( INCX.EQ.1 )THEN
2668             DO 20, J = 1, N
2669                IF( X( J ).NE.ZERO )THEN
2670                   TEMP = ALPHA*X( J )
2671                   DO 10, I = 1, J
2672                      A( I, J ) = A( I, J ) + X( I )*TEMP
2673    10             CONTINUE
2674                END IF
2675    20       CONTINUE
2676          ELSE
2677             JX = KX
2678             DO 40, J = 1, N
2679                IF( X( JX ).NE.ZERO )THEN
2680                   TEMP = ALPHA*X( JX )
2681                   IX   = KX
2682                   DO 30, I = 1, J
2683                      A( I, J ) = A( I, J ) + X( IX )*TEMP
2684                      IX        = IX        + INCX
2685    30             CONTINUE
2686                END IF
2687                JX = JX + INCX
2688    40       CONTINUE
2689          END IF
2690       ELSE
2691 *
2692 *        Form  A  when A is stored in lower triangle.
2693 *
2694          IF( INCX.EQ.1 )THEN
2695             DO 60, J = 1, N
2696                IF( X( J ).NE.ZERO )THEN
2697                   TEMP = ALPHA*X( J )
2698                   DO 50, I = J, N
2699                      A( I, J ) = A( I, J ) + X( I )*TEMP
2700    50             CONTINUE
2701                END IF
2702    60       CONTINUE
2703          ELSE
2704             JX = KX
2705             DO 80, J = 1, N
2706                IF( X( JX ).NE.ZERO )THEN
2707                   TEMP = ALPHA*X( JX )
2708                   IX   = JX
2709                   DO 70, I = J, N
2710                      A( I, J ) = A( I, J ) + X( IX )*TEMP
2711                      IX        = IX        + INCX
2712    70             CONTINUE
2713                END IF
2714                JX = JX + INCX
2715    80       CONTINUE
2716          END IF
2717       END IF
2718 *
2719       RETURN
2720 *
2721 *     End of DSYR  .
2722 *
2723       END
2724       SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
2725      $                   BETA, C, LDC )
2726 *     .. Scalar Arguments ..
2727       CHARACTER*1        UPLO, TRANS
2728       INTEGER            N, K, LDA, LDC
2729       DOUBLE PRECISION   ALPHA, BETA
2730 *     .. Array Arguments ..
2731       DOUBLE PRECISION   A( LDA, * ), C( LDC, * )
2732 *     ..
2733 *
2734 *  Purpose
2735 *  =======
2736 *
2737 *  DSYRK  performs one of the symmetric rank k operations
2738 *
2739 *     C := alpha*A*A' + beta*C,
2740 *
2741 *  or
2742 *
2743 *     C := alpha*A'*A + beta*C,
2744 *
2745 *  where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
2746 *  and  A  is an  n by k  matrix in the first case and a  k by n  matrix
2747 *  in the second case.
2748 *
2749 *  Parameters
2750 *  ==========
2751 *
2752 *  UPLO   - CHARACTER*1.
2753 *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
2754 *           triangular  part  of the  array  C  is to be  referenced  as
2755 *           follows:
2756 *
2757 *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
2758 *                                  is to be referenced.
2759 *
2760 *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
2761 *                                  is to be referenced.
2762 *
2763 *           Unchanged on exit.
2764 *
2765 *  TRANS  - CHARACTER*1.
2766 *           On entry,  TRANS  specifies the operation to be performed as
2767 *           follows:
2768 *
2769 *              TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C.
2770 *
2771 *              TRANS = 'T' or 't'   C := alpha*A'*A + beta*C.
2772 *
2773 *              TRANS = 'C' or 'c'   C := alpha*A'*A + beta*C.
2774 *
2775 *           Unchanged on exit.
2776 *
2777 *  N      - INTEGER.
2778 *           On entry,  N specifies the order of the matrix C.  N must be
2779 *           at least zero.
2780 *           Unchanged on exit.
2781 *
2782 *  K      - INTEGER.
2783 *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
2784 *           of  columns   of  the   matrix   A,   and  on   entry   with
2785 *           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
2786 *           of rows of the matrix  A.  K must be at least zero.
2787 *           Unchanged on exit.
2788 *
2789 *  ALPHA  - DOUBLE PRECISION.
2790 *           On entry, ALPHA specifies the scalar alpha.
2791 *           Unchanged on exit.
2792 *
2793 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
2794 *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
2795 *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
2796 *           part of the array  A  must contain the matrix  A,  otherwise
2797 *           the leading  k by n  part of the array  A  must contain  the
2798 *           matrix A.
2799 *           Unchanged on exit.
2800 *
2801 *  LDA    - INTEGER.
2802 *           On entry, LDA specifies the first dimension of A as declared
2803 *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
2804 *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
2805 *           be at least  max( 1, k ).
2806 *           Unchanged on exit.
2807 *
2808 *  BETA   - DOUBLE PRECISION.
2809 *           On entry, BETA specifies the scalar beta.
2810 *           Unchanged on exit.
2811 *
2812 *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
2813 *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
2814 *           upper triangular part of the array C must contain the upper
2815 *           triangular part  of the  symmetric matrix  and the strictly
2816 *           lower triangular part of C is not referenced.  On exit, the
2817 *           upper triangular part of the array  C is overwritten by the
2818 *           upper triangular part of the updated matrix.
2819 *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
2820 *           lower triangular part of the array C must contain the lower
2821 *           triangular part  of the  symmetric matrix  and the strictly
2822 *           upper triangular part of C is not referenced.  On exit, the
2823 *           lower triangular part of the array  C is overwritten by the
2824 *           lower triangular part of the updated matrix.
2825 *
2826 *  LDC    - INTEGER.
2827 *           On entry, LDC specifies the first dimension of C as declared
2828 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
2829 *           max( 1, n ).
2830 *           Unchanged on exit.
2831 *
2832 *
2833 *  Level 3 Blas routine.
2834 *
2835 *  -- Written on 8-February-1989.
2836 *     Jack Dongarra, Argonne National Laboratory.
2837 *     Iain Duff, AERE Harwell.
2838 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
2839 *     Sven Hammarling, Numerical Algorithms Group Ltd.
2840 *
2841 *
2842 *     .. External Functions ..
2843       LOGICAL            LSAME
2844       EXTERNAL           LSAME
2845 *     .. External Subroutines ..
2846       EXTERNAL           XERBLA
2847 *     .. Intrinsic Functions ..
2848       INTRINSIC          MAX
2849 *     .. Local Scalars ..
2850       LOGICAL            UPPER
2851       INTEGER            I, INFO, J, L, NROWA
2852       DOUBLE PRECISION   TEMP
2853 *     .. Parameters ..
2854       DOUBLE PRECISION   ONE ,         ZERO
2855       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
2856 *     ..
2857 *     .. Executable Statements ..
2858 *
2859 *     Test the input parameters.
2860 *
2861       IF( LSAME( TRANS, 'N' ) )THEN
2862          NROWA = N
2863       ELSE
2864          NROWA = K
2865       END IF
2866       UPPER = LSAME( UPLO, 'U' )
2867 *
2868       INFO = 0
2869       IF(      ( .NOT.UPPER               ).AND.
2870      $         ( .NOT.LSAME( UPLO , 'L' ) )      )THEN
2871          INFO = 1
2872       ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
2873      $         ( .NOT.LSAME( TRANS, 'T' ) ).AND.
2874      $         ( .NOT.LSAME( TRANS, 'C' ) )      )THEN
2875          INFO = 2
2876       ELSE IF( N  .LT.0               )THEN
2877          INFO = 3
2878       ELSE IF( K  .LT.0               )THEN
2879          INFO = 4
2880       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
2881          INFO = 7
2882       ELSE IF( LDC.LT.MAX( 1, N     ) )THEN
2883          INFO = 10
2884       END IF
2885       IF( INFO.NE.0 )THEN
2886          CALL XERBLA( 'DSYRK ', INFO )
2887          RETURN
2888       END IF
2889 *
2890 *     Quick return if possible.
2891 *
2892       IF( ( N.EQ.0 ).OR.
2893      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
2894      $   RETURN
2895 *
2896 *     And when  alpha.eq.zero.
2897 *
2898       IF( ALPHA.EQ.ZERO )THEN
2899          IF( UPPER )THEN
2900             IF( BETA.EQ.ZERO )THEN
2901                DO 20, J = 1, N
2902                   DO 10, I = 1, J
2903                      C( I, J ) = ZERO
2904    10             CONTINUE
2905    20          CONTINUE
2906             ELSE
2907                DO 40, J = 1, N
2908                   DO 30, I = 1, J
2909                      C( I, J ) = BETA*C( I, J )
2910    30             CONTINUE
2911    40          CONTINUE
2912             END IF
2913          ELSE
2914             IF( BETA.EQ.ZERO )THEN
2915                DO 60, J = 1, N
2916                   DO 50, I = J, N
2917                      C( I, J ) = ZERO
2918    50             CONTINUE
2919    60          CONTINUE
2920             ELSE
2921                DO 80, J = 1, N
2922                   DO 70, I = J, N
2923                      C( I, J ) = BETA*C( I, J )
2924    70             CONTINUE
2925    80          CONTINUE
2926             END IF
2927          END IF
2928          RETURN
2929       END IF
2930 *
2931 *     Start the operations.
2932 *
2933       IF( LSAME( TRANS, 'N' ) )THEN
2934 *
2935 *        Form  C := alpha*A*A' + beta*C.
2936 *
2937          IF( UPPER )THEN
2938             DO 130, J = 1, N
2939                IF( BETA.EQ.ZERO )THEN
2940                   DO 90, I = 1, J
2941                      C( I, J ) = ZERO
2942    90             CONTINUE
2943                ELSE IF( BETA.NE.ONE )THEN
2944                   DO 100, I = 1, J
2945                      C( I, J ) = BETA*C( I, J )
2946   100             CONTINUE
2947                END IF
2948                DO 120, L = 1, K
2949                   IF( A( J, L ).NE.ZERO )THEN
2950                      TEMP = ALPHA*A( J, L )
2951                      DO 110, I = 1, J
2952                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
2953   110                CONTINUE
2954                   END IF
2955   120          CONTINUE
2956   130       CONTINUE
2957          ELSE
2958             DO 180, J = 1, N
2959                IF( BETA.EQ.ZERO )THEN
2960                   DO 140, I = J, N
2961                      C( I, J ) = ZERO
2962   140             CONTINUE
2963                ELSE IF( BETA.NE.ONE )THEN
2964                   DO 150, I = J, N
2965                      C( I, J ) = BETA*C( I, J )
2966   150             CONTINUE
2967                END IF
2968                DO 170, L = 1, K
2969                   IF( A( J, L ).NE.ZERO )THEN
2970                      TEMP      = ALPHA*A( J, L )
2971                      DO 160, I = J, N
2972                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
2973   160                CONTINUE
2974                   END IF
2975   170          CONTINUE
2976   180       CONTINUE
2977          END IF
2978       ELSE
2979 *
2980 *        Form  C := alpha*A'*A + beta*C.
2981 *
2982          IF( UPPER )THEN
2983             DO 210, J = 1, N
2984                DO 200, I = 1, J
2985                   TEMP = ZERO
2986                   DO 190, L = 1, K
2987                      TEMP = TEMP + A( L, I )*A( L, J )
2988   190             CONTINUE
2989                   IF( BETA.EQ.ZERO )THEN
2990                      C( I, J ) = ALPHA*TEMP
2991                   ELSE
2992                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
2993                   END IF
2994   200          CONTINUE
2995   210       CONTINUE
2996          ELSE
2997             DO 240, J = 1, N
2998                DO 230, I = J, N
2999                   TEMP = ZERO
3000                   DO 220, L = 1, K
3001                      TEMP = TEMP + A( L, I )*A( L, J )
3002   220             CONTINUE
3003                   IF( BETA.EQ.ZERO )THEN
3004                      C( I, J ) = ALPHA*TEMP
3005                   ELSE
3006                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
3007                   END IF
3008   230          CONTINUE
3009   240       CONTINUE
3010          END IF
3011       END IF
3012 *
3013       RETURN
3014 *
3015 *     End of DSYRK .
3016 *
3017       END
3018       SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
3019      $                   BETA, C, LDC )
3020 *     .. Scalar Arguments ..
3021       CHARACTER*1        TRANSA, TRANSB
3022       INTEGER            M, N, K, LDA, LDB, LDC
3023       DOUBLE PRECISION   ALPHA, BETA
3024 *     .. Array Arguments ..
3025       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
3026 *     ..
3027 *
3028 *  Purpose
3029 *  =======
3030 *
3031 *  DGEMM  performs one of the matrix-matrix operations
3032 *
3033 *     C := alpha*op( A )*op( B ) + beta*C,
3034 *
3035 *  where  op( X ) is one of
3036 *
3037 *     op( X ) = X   or   op( X ) = X',
3038 *
3039 *  alpha and beta are scalars, and A, B and C are matrices, with op( A )
3040 *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
3041 *
3042 *  Parameters
3043 *  ==========
3044 *
3045 *  TRANSA - CHARACTER*1.
3046 *           On entry, TRANSA specifies the form of op( A ) to be used in
3047 *           the matrix multiplication as follows:
3048 *
3049 *              TRANSA = 'N' or 'n',  op( A ) = A.
3050 *
3051 *              TRANSA = 'T' or 't',  op( A ) = A'.
3052 *
3053 *              TRANSA = 'C' or 'c',  op( A ) = A'.
3054 *
3055 *           Unchanged on exit.
3056 *
3057 *  TRANSB - CHARACTER*1.
3058 *           On entry, TRANSB specifies the form of op( B ) to be used in
3059 *           the matrix multiplication as follows:
3060 *
3061 *              TRANSB = 'N' or 'n',  op( B ) = B.
3062 *
3063 *              TRANSB = 'T' or 't',  op( B ) = B'.
3064 *
3065 *              TRANSB = 'C' or 'c',  op( B ) = B'.
3066 *
3067 *           Unchanged on exit.
3068 *
3069 *  M      - INTEGER.
3070 *           On entry,  M  specifies  the number  of rows  of the  matrix
3071 *           op( A )  and of the  matrix  C.  M  must  be at least  zero.
3072 *           Unchanged on exit.
3073 *
3074 *  N      - INTEGER.
3075 *           On entry,  N  specifies the number  of columns of the matrix
3076 *           op( B ) and the number of columns of the matrix C. N must be
3077 *           at least zero.
3078 *           Unchanged on exit.
3079 *
3080 *  K      - INTEGER.
3081 *           On entry,  K  specifies  the number of columns of the matrix
3082 *           op( A ) and the number of rows of the matrix op( B ). K must
3083 *           be at least  zero.
3084 *           Unchanged on exit.
3085 *
3086 *  ALPHA  - DOUBLE PRECISION.
3087 *           On entry, ALPHA specifies the scalar alpha.
3088 *           Unchanged on exit.
3089 *
3090 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
3091 *           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
3092 *           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
3093 *           part of the array  A  must contain the matrix  A,  otherwise
3094 *           the leading  k by m  part of the array  A  must contain  the
3095 *           matrix A.
3096 *           Unchanged on exit.
3097 *
3098 *  LDA    - INTEGER.
3099 *           On entry, LDA specifies the first dimension of A as declared
3100 *           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
3101 *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
3102 *           least  max( 1, k ).
3103 *           Unchanged on exit.
3104 *
3105 *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
3106 *           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
3107 *           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
3108 *           part of the array  B  must contain the matrix  B,  otherwise
3109 *           the leading  n by k  part of the array  B  must contain  the
3110 *           matrix B.
3111 *           Unchanged on exit.
3112 *
3113 *  LDB    - INTEGER.
3114 *           On entry, LDB specifies the first dimension of B as declared
3115 *           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
3116 *           LDB must be at least  max( 1, k ), otherwise  LDB must be at
3117 *           least  max( 1, n ).
3118 *           Unchanged on exit.
3119 *
3120 *  BETA   - DOUBLE PRECISION.
3121 *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
3122 *           supplied as zero then C need not be set on input.
3123 *           Unchanged on exit.
3124 *
3125 *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
3126 *           Before entry, the leading  m by n  part of the array  C must
3127 *           contain the matrix  C,  except when  beta  is zero, in which
3128 *           case C need not be set on entry.
3129 *           On exit, the array  C  is overwritten by the  m by n  matrix
3130 *           ( alpha*op( A )*op( B ) + beta*C ).
3131 *
3132 *  LDC    - INTEGER.
3133 *           On entry, LDC specifies the first dimension of C as declared
3134 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
3135 *           max( 1, m ).
3136 *           Unchanged on exit.
3137 *
3138 *
3139 *  Level 3 Blas routine.
3140 *
3141 *  -- Written on 8-February-1989.
3142 *     Jack Dongarra, Argonne National Laboratory.
3143 *     Iain Duff, AERE Harwell.
3144 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
3145 *     Sven Hammarling, Numerical Algorithms Group Ltd.
3146 *
3147 *
3148 *     .. External Functions ..
3149       LOGICAL            LSAME
3150       EXTERNAL           LSAME
3151 *     .. External Subroutines ..
3152       EXTERNAL           XERBLA
3153 *     .. Intrinsic Functions ..
3154       INTRINSIC          MAX
3155 *     .. Local Scalars ..
3156       LOGICAL            NOTA, NOTB
3157       INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
3158       DOUBLE PRECISION   TEMP
3159 *     .. Parameters ..
3160       DOUBLE PRECISION   ONE         , ZERO
3161       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
3162 *     ..
3163 *     .. Executable Statements ..
3164 *
3165 *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
3166 *     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
3167 *     and  columns of  A  and the  number of  rows  of  B  respectively.
3168 *
3169       NOTA  = LSAME( TRANSA, 'N' )
3170       NOTB  = LSAME( TRANSB, 'N' )
3171       IF( NOTA )THEN
3172          NROWA = M
3173          NCOLA = K
3174       ELSE
3175          NROWA = K
3176          NCOLA = M
3177       END IF
3178       IF( NOTB )THEN
3179          NROWB = K
3180       ELSE
3181          NROWB = N
3182       END IF
3183 *
3184 *     Test the input parameters.
3185 *
3186       INFO = 0
3187       IF(      ( .NOT.NOTA                 ).AND.
3188      $         ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
3189      $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
3190          INFO = 1
3191       ELSE IF( ( .NOT.NOTB                 ).AND.
3192      $         ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
3193      $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
3194          INFO = 2
3195       ELSE IF( M  .LT.0               )THEN
3196          INFO = 3
3197       ELSE IF( N  .LT.0               )THEN
3198          INFO = 4
3199       ELSE IF( K  .LT.0               )THEN
3200          INFO = 5
3201       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
3202          INFO = 8
3203       ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
3204          INFO = 10
3205       ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
3206          INFO = 13
3207       END IF
3208       IF( INFO.NE.0 )THEN
3209          CALL XERBLA( 'DGEMM ', INFO )
3210          RETURN
3211       END IF
3212 *
3213 *     Quick return if possible.
3214 *
3215       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
3216      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
3217      $   RETURN
3218 *
3219 *     And if  alpha.eq.zero.
3220 *
3221       IF( ALPHA.EQ.ZERO )THEN
3222          IF( BETA.EQ.ZERO )THEN
3223             DO 20, J = 1, N
3224                DO 10, I = 1, M
3225                   C( I, J ) = ZERO
3226    10          CONTINUE
3227    20       CONTINUE
3228          ELSE
3229             DO 40, J = 1, N
3230                DO 30, I = 1, M
3231                   C( I, J ) = BETA*C( I, J )
3232    30          CONTINUE
3233    40       CONTINUE
3234          END IF
3235          RETURN
3236       END IF
3237 *
3238 *     Start the operations.
3239 *
3240       IF( NOTB )THEN
3241          IF( NOTA )THEN
3242 *
3243 *           Form  C := alpha*A*B + beta*C.
3244 *
3245             DO 90, J = 1, N
3246                IF( BETA.EQ.ZERO )THEN
3247                   DO 50, I = 1, M
3248                      C( I, J ) = ZERO
3249    50             CONTINUE
3250                ELSE IF( BETA.NE.ONE )THEN
3251                   DO 60, I = 1, M
3252                      C( I, J ) = BETA*C( I, J )
3253    60             CONTINUE
3254                END IF
3255                DO 80, L = 1, K
3256                   IF( B( L, J ).NE.ZERO )THEN
3257                      TEMP = ALPHA*B( L, J )
3258                      DO 70, I = 1, M
3259                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
3260    70                CONTINUE
3261                   END IF
3262    80          CONTINUE
3263    90       CONTINUE
3264          ELSE
3265 *
3266 *           Form  C := alpha*A'*B + beta*C
3267 *
3268             DO 120, J = 1, N
3269                DO 110, I = 1, M
3270                   TEMP = ZERO
3271                   DO 100, L = 1, K
3272                      TEMP = TEMP + A( L, I )*B( L, J )
3273   100             CONTINUE
3274                   IF( BETA.EQ.ZERO )THEN
3275                      C( I, J ) = ALPHA*TEMP
3276                   ELSE
3277                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
3278                   END IF
3279   110          CONTINUE
3280   120       CONTINUE
3281          END IF
3282       ELSE
3283          IF( NOTA )THEN
3284 *
3285 *           Form  C := alpha*A*B' + beta*C
3286 *
3287             DO 170, J = 1, N
3288                IF( BETA.EQ.ZERO )THEN
3289                   DO 130, I = 1, M
3290                      C( I, J ) = ZERO
3291   130             CONTINUE
3292                ELSE IF( BETA.NE.ONE )THEN
3293                   DO 140, I = 1, M
3294                      C( I, J ) = BETA*C( I, J )
3295   140             CONTINUE
3296                END IF
3297                DO 160, L = 1, K
3298                   IF( B( J, L ).NE.ZERO )THEN
3299                      TEMP = ALPHA*B( J, L )
3300                      DO 150, I = 1, M
3301                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
3302   150                CONTINUE
3303                   END IF
3304   160          CONTINUE
3305   170       CONTINUE
3306          ELSE
3307 *
3308 *           Form  C := alpha*A'*B' + beta*C
3309 *
3310             DO 200, J = 1, N
3311                DO 190, I = 1, M
3312                   TEMP = ZERO
3313                   DO 180, L = 1, K
3314                      TEMP = TEMP + A( L, I )*B( J, L )
3315   180             CONTINUE
3316                   IF( BETA.EQ.ZERO )THEN
3317                      C( I, J ) = ALPHA*TEMP
3318                   ELSE
3319                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
3320                   END IF
3321   190          CONTINUE
3322   200       CONTINUE
3323          END IF
3324       END IF
3325 *
3326       RETURN
3327 *
3328 *     End of DGEMM .
3329 *
3330       END
3331       SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
3332      $                   BETA, Y, INCY )
3333 *     .. Scalar Arguments ..
3334       DOUBLE PRECISION   ALPHA, BETA
3335       INTEGER            INCX, INCY, LDA, M, N
3336       CHARACTER*1        TRANS
3337 *     .. Array Arguments ..
3338       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
3339 *     ..
3340 *
3341 *  Purpose
3342 *  =======
3343 *
3344 *  DGEMV  performs one of the matrix-vector operations
3345 *
3346 *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
3347 *
3348 *  where alpha and beta are scalars, x and y are vectors and A is an
3349 *  m by n matrix.
3350 *
3351 *  Parameters
3352 *  ==========
3353 *
3354 *  TRANS  - CHARACTER*1.
3355 *           On entry, TRANS specifies the operation to be performed as
3356 *           follows:
3357 *
3358 *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
3359 *
3360 *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
3361 *
3362 *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
3363 *
3364 *           Unchanged on exit.
3365 *
3366 *  M      - INTEGER.
3367 *           On entry, M specifies the number of rows of the matrix A.
3368 *           M must be at least zero.
3369 *           Unchanged on exit.
3370 *
3371 *  N      - INTEGER.
3372 *           On entry, N specifies the number of columns of the matrix A.
3373 *           N must be at least zero.
3374 *           Unchanged on exit.
3375 *
3376 *  ALPHA  - DOUBLE PRECISION.
3377 *           On entry, ALPHA specifies the scalar alpha.
3378 *           Unchanged on exit.
3379 *
3380 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
3381 *           Before entry, the leading m by n part of the array A must
3382 *           contain the matrix of coefficients.
3383 *           Unchanged on exit.
3384 *
3385 *  LDA    - INTEGER.
3386 *           On entry, LDA specifies the first dimension of A as declared
3387 *           in the calling (sub) program. LDA must be at least
3388 *           max( 1, m ).
3389 *           Unchanged on exit.
3390 *
3391 *  X      - DOUBLE PRECISION array of DIMENSION at least
3392 *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
3393 *           and at least
3394 *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
3395 *           Before entry, the incremented array X must contain the
3396 *           vector x.
3397 *           Unchanged on exit.
3398 *
3399 *  INCX   - INTEGER.
3400 *           On entry, INCX specifies the increment for the elements of
3401 *           X. INCX must not be zero.
3402 *           Unchanged on exit.
3403 *
3404 *  BETA   - DOUBLE PRECISION.
3405 *           On entry, BETA specifies the scalar beta. When BETA is
3406 *           supplied as zero then Y need not be set on input.
3407 *           Unchanged on exit.
3408 *
3409 *  Y      - DOUBLE PRECISION array of DIMENSION at least
3410 *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
3411 *           and at least
3412 *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
3413 *           Before entry with BETA non-zero, the incremented array Y
3414 *           must contain the vector y. On exit, Y is overwritten by the
3415 *           updated vector y.
3416 *
3417 *  INCY   - INTEGER.
3418 *           On entry, INCY specifies the increment for the elements of
3419 *           Y. INCY must not be zero.
3420 *           Unchanged on exit.
3421 *
3422 *
3423 *  Level 2 Blas routine.
3424 *
3425 *  -- Written on 22-October-1986.
3426 *     Jack Dongarra, Argonne National Lab.
3427 *     Jeremy Du Croz, Nag Central Office.
3428 *     Sven Hammarling, Nag Central Office.
3429 *     Richard Hanson, Sandia National Labs.
3430 *
3431 *
3432 *     .. Parameters ..
3433       DOUBLE PRECISION   ONE         , ZERO
3434       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
3435 *     .. Local Scalars ..
3436       DOUBLE PRECISION   TEMP
3437       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
3438 *     .. External Functions ..
3439       LOGICAL            LSAME
3440       EXTERNAL           LSAME
3441 *     .. External Subroutines ..
3442       EXTERNAL           XERBLA
3443 *     .. Intrinsic Functions ..
3444       INTRINSIC          MAX
3445 *     ..
3446 *     .. Executable Statements ..
3447 *
3448 *     Test the input parameters.
3449 *
3450       INFO = 0
3451       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
3452      $         .NOT.LSAME( TRANS, 'T' ).AND.
3453      $         .NOT.LSAME( TRANS, 'C' )      )THEN
3454          INFO = 1
3455       ELSE IF( M.LT.0 )THEN
3456          INFO = 2
3457       ELSE IF( N.LT.0 )THEN
3458          INFO = 3
3459       ELSE IF( LDA.LT.MAX( 1, M ) )THEN
3460          INFO = 6
3461       ELSE IF( INCX.EQ.0 )THEN
3462          INFO = 8
3463       ELSE IF( INCY.EQ.0 )THEN
3464          INFO = 11
3465       END IF
3466       IF( INFO.NE.0 )THEN
3467          CALL XERBLA( 'DGEMV ', INFO )
3468          RETURN
3469       END IF
3470 *
3471 *     Quick return if possible.
3472 *
3473       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
3474      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
3475      $   RETURN
3476 *
3477 *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
3478 *     up the start points in  X  and  Y.
3479 *
3480       IF( LSAME( TRANS, 'N' ) )THEN
3481          LENX = N
3482          LENY = M
3483       ELSE
3484          LENX = M
3485          LENY = N
3486       END IF
3487       IF( INCX.GT.0 )THEN
3488          KX = 1
3489       ELSE
3490          KX = 1 - ( LENX - 1 )*INCX
3491       END IF
3492       IF( INCY.GT.0 )THEN
3493          KY = 1
3494       ELSE
3495          KY = 1 - ( LENY - 1 )*INCY
3496       END IF
3497 *
3498 *     Start the operations. In this version the elements of A are
3499 *     accessed sequentially with one pass through A.
3500 *
3501 *     First form  y := beta*y.
3502 *
3503       IF( BETA.NE.ONE )THEN
3504          IF( INCY.EQ.1 )THEN
3505             IF( BETA.EQ.ZERO )THEN
3506                DO 10, I = 1, LENY
3507                   Y( I ) = ZERO
3508    10          CONTINUE
3509             ELSE
3510                DO 20, I = 1, LENY
3511                   Y( I ) = BETA*Y( I )
3512    20          CONTINUE
3513             END IF
3514          ELSE
3515             IY = KY
3516             IF( BETA.EQ.ZERO )THEN
3517                DO 30, I = 1, LENY
3518                   Y( IY ) = ZERO
3519                   IY      = IY   + INCY
3520    30          CONTINUE
3521             ELSE
3522                DO 40, I = 1, LENY
3523                   Y( IY ) = BETA*Y( IY )
3524                   IY      = IY           + INCY
3525    40          CONTINUE
3526             END IF
3527          END IF
3528       END IF
3529       IF( ALPHA.EQ.ZERO )
3530      $   RETURN
3531       IF( LSAME( TRANS, 'N' ) )THEN
3532 *
3533 *        Form  y := alpha*A*x + y.
3534 *
3535          JX = KX
3536          IF( INCY.EQ.1 )THEN
3537             DO 60, J = 1, N
3538                IF( X( JX ).NE.ZERO )THEN
3539                   TEMP = ALPHA*X( JX )
3540                   DO 50, I = 1, M
3541                      Y( I ) = Y( I ) + TEMP*A( I, J )
3542    50             CONTINUE
3543                END IF
3544                JX = JX + INCX
3545    60       CONTINUE
3546          ELSE
3547             DO 80, J = 1, N
3548                IF( X( JX ).NE.ZERO )THEN
3549                   TEMP = ALPHA*X( JX )
3550                   IY   = KY
3551                   DO 70, I = 1, M
3552                      Y( IY ) = Y( IY ) + TEMP*A( I, J )
3553                      IY      = IY      + INCY
3554    70             CONTINUE
3555                END IF
3556                JX = JX + INCX
3557    80       CONTINUE
3558          END IF
3559       ELSE
3560 *
3561 *        Form  y := alpha*A'*x + y.
3562 *
3563          JY = KY
3564          IF( INCX.EQ.1 )THEN
3565             DO 100, J = 1, N
3566                TEMP = ZERO
3567                DO 90, I = 1, M
3568                   TEMP = TEMP + A( I, J )*X( I )
3569    90          CONTINUE
3570                Y( JY ) = Y( JY ) + ALPHA*TEMP
3571                JY      = JY      + INCY
3572   100       CONTINUE
3573          ELSE
3574             DO 120, J = 1, N
3575                TEMP = ZERO
3576                IX   = KX
3577                DO 110, I = 1, M
3578                   TEMP = TEMP + A( I, J )*X( IX )
3579                   IX   = IX   + INCX
3580   110          CONTINUE
3581                Y( JY ) = Y( JY ) + ALPHA*TEMP
3582                JY      = JY      + INCY
3583   120       CONTINUE
3584          END IF
3585       END IF
3586 *
3587       RETURN
3588 *
3589 *     End of DGEMV .
3590 *
3591       END
3592
3593       subroutine perma(compo)
3594       include "calcium.hf"
3595       double precision a(5,24),b(24),q(6)
3596       dimension t(24),tn(24)
3597       dimension puissn(24),puiss(24)
3598       dimension text(6),rflu(6),textn(6)
3599       dimension tpar(6),rpar(6),tparn(6)
3600       integer compo
3601
3602       ldab=5
3603       ldb=24
3604       nrhs=1
3605       npt=0
3606       r=1.
3607       ro=1.
3608       rext=0.5
3609       te=10.
3610       do i=1,24
3611       tn(i)=100.
3612       t(i)=100.
3613       puiss(i)=100.
3614       enddo
3615       do j=1,6
3616       q(j)=50.
3617       enddo
3618 c
3619 c   construction de la matrice (laplacien)
3620 c
3621       do i = 1, 5
3622         do j = 1, 24
3623           a(i,j) = 0.
3624         enddo
3625       enddo
3626
3627       do i=2,3
3628       do j=2,5
3629       npt=i+(j-1)*4
3630       a(1,npt)=4./r
3631       a(2,npt)=-1./r
3632       a(5,npt)=-1./r
3633       enddo
3634       enddo
3635       do i=2,3
3636       npt=i
3637       a(1,npt)=3./r
3638       a(2,npt)=-1./r
3639       a(5,npt)=-1./r
3640       npt=i+20
3641       a(1,npt)=3./r
3642       a(2,npt)=-1./r
3643       a(5,npt)=0.
3644       enddo
3645       do j=2,5
3646       npt=1+(j-1)*4
3647       a(1,npt)=3./r
3648       a(2,npt)=-1./r
3649       a(5,npt)=-1./r
3650       npt=4+(j-1)*4
3651       a(1,npt)=3./r+1./(r/2.+rext)
3652       a(2,npt)=0.
3653       a(5,npt)=-1./r
3654       enddo
3655       i=1
3656       a(1,i)=2./r
3657       a(2,i)=-1./r
3658       a(5,i)=-1./r
3659       i=21
3660       a(1,i)=2./r
3661       a(2,i)=-1./r
3662       a(5,i)=-1./r
3663       i=24
3664       a(1,i)=2./r+1./(r/2.+rext)
3665       a(2,i)=-1./r
3666       a(5,i)=-1./r
3667       i=4
3668       a(1,i)=2./r+1./(r/2.+rext)
3669       a(2,i)=0.
3670       a(5,i)=-1./r
3671
3672 c
3673 c  factorisation de la matrice
3674 c
3675       n=24
3676       kd=4
3677
3678       call dPBTRF( 'L' , N, KD, A, LDAB, INFO )
3679
3680       iter=0
3681
3682       do i=1,6
3683       tpar(i)=t(4*i)
3684       rpar(i)=r/2.
3685       enddo
3686
3687       do i=1,24
3688       tn(i)=0.
3689       puissn(i)=0.
3690       enddo
3691       do i=1,6
3692       tparn(i)=0.
3693       textn(i)=0.
3694       rflu(i)=0.
3695       enddo
3696 c
3697 c  initialisation de la temperature a iter=0
3698 c
3699         CALL cpeRE(compo,CP_ITERATION, ti, iter, 'tempi', 24,
3700      &            t , info)
3701             IF( info.EQ.CPSTOP )GO TO 9000
3702 c
3703 c  initialisation de la temperature de bord a iter=0.
3704 c
3705         CALL cpeRE(compo,CP_ITERATION, ti, iter, 'tpi', 6,
3706      &            tpar , info)
3707             IF( info.EQ.CPSTOP )GO TO 9000
3708 c
3709 c  boucle iterative  infinie
3710 c
3711       do while( iter .lt. 500)
3712 c
3713 c  lecture de la puissance combustible a iter
3714 c
3715         CALL cplRE(compo,CP_ITERATION,ti, tf, iter, 'puissi', 24,
3716      &      nval, puiss , info)
3717             IF( info.EQ.CPSTOP )GO TO 9000
3718 c
3719 c  lecture de la temperature exterieure a iter
3720 c
3721         CALL cplRE(compo,CP_ITERATION,ti, tf, iter, 'tfi', 6,
3722      &      nval, text , info)
3723             IF( info.EQ.CPSTOP )GO TO 9000
3724 c
3725 c  calcul du second membre
3726 c
3727       do npt=1,24
3728       b(npt)=puiss(npt)
3729       enddo
3730
3731       do j=1,6
3732       npt=j*4
3733       b(npt)=b(npt)+text(j)/(r/2+rflu(j))
3734       enddo
3735
3736 c
3737 c  resolution du systeme lineaire
3738 c
3739       call dPBTRs( 'L' , N,kd, nrhs,A,LDAB,b,ldb,INFO )
3740
3741       do i=1,24
3742       t(i)=b(i)
3743       enddo
3744
3745       do i=1,6
3746       tpar(i)=t(4*i)
3747       enddo
3748
3749       iter=iter+1
3750 c
3751 c  ecriture de la temperature a iter+1
3752 c
3753         CALL cpeRE(compo,CP_ITERATION, ti, iter, 'tempi', 24,
3754      &            t , info)
3755             IF( info.NE. CPOK  )GO TO 9000
3756 c
3757 c  ecriture de la temperature de paroi a iter+1
3758 c
3759         CALL cpeRE(compo,CP_ITERATION, ti, iter, 'tpi', 6,
3760      &            tpar , info)
3761             IF( info.NE. CPOK  )GO TO 9000
3762 c
3763 c  calcul de convergence
3764 c
3765       conv1=0.
3766       do i=1,npt
3767       conv1=conv1+(puiss(i)-puissn(i))**2
3768       conv1=conv1+(t(i)-tn(i))**2
3769       enddo
3770       do i=1,6
3771       conv1=conv1+(tpar(i)-tparn(i))**2
3772       conv1=conv1+(text(i)-textn(i))**2
3773       enddo
3774
3775       iconv=0
3776       if(conv1.lt.1.e-3) iconv=1
3777 c
3778 c   ecriture du flag de convergence iconv
3779 c
3780       write(6,*)"SOLIDE:",iter,iconv
3781       call flush(6)
3782         CALL cpeEN(compo,CP_ITERATION,tf, iter  , 'iconv', 1,
3783      &     iconv  , info)
3784             IF( info.NE. CPOK  )GO TO 9000
3785
3786       if(iconv.eq.1)go to 9000
3787
3788       do i=1,24
3789       tn(i)=b(i)
3790       puissn(i)=puiss(i)
3791       enddo
3792       do i=1,6
3793       tparn(i)=tpar(i)
3794       textn(i)=text(i)
3795       enddo
3796
3797       enddo
3798
3799 9000  continue
3800             CALL cpfin(compo,CP_ARRET, info)
3801       end