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