1 subroutine transit(compo,dt0)
3 double precision a(5,24),b(24),tn(24)
5 dimension maille(2),t(24),puiss(24)
6 dimension text(6),rflu(6)
7 dimension tpar(6),rpar(6)
31 c construction de la matrice (laplacien)
63 a(1,npt)=3./r+1./(r/2.+rext)+ro/dt
76 a(1,i)=2./r+1./(r/2.+rext)+ro/dt
80 a(1,i)=2./r+1./(r/2.+rext)+ro/dt
84 c factorisation de la matrice
89 call dPBTRF( 'L' , N, KD, A, LDAB, INFO )
101 c initialisation de la temperature a t=0
103 CALL cpeRE(compo,CP_TEMPS, ti, npas, 'temperature', 24,
105 IF( info.NE. CPOK )GO TO 9000
107 c initialisation de la temperature de bord a t=0.
109 CALL cpeRE(compo,CP_TEMPS, ti, npas, 'tparoi', 6,
111 IF( info.NE. CPOK )GO TO 9000
113 c initialisation de la resistance thermique de bord a t=0.
115 CALL cpeRE(compo,CP_TEMPS, ti, npas, 'rparoi', 6,
117 IF( info.NE. CPOK )GO TO 9000
119 c boucle temporelle infinie
122 c do while( ti.lT.100. )
125 c lecture de la puissance combustible entre ti et tf
127 CALL cplRE(compo,CP_TEMPS,ti, tf, npas, 'puissa', 24,
128 & nval, puiss , info)
129 IF( info.NE. CPOK )GO TO 9000
131 c lecture de la temperature exterieure entre ti et tf
133 CALL cplRE(compo,CP_TEMPS,ti, tf, npas, 'text', 6,
135 IF( info.NE. CPOK )GO TO 9000
137 c lecture de la resistance exterieure entre ti et tf
139 CALL cplRE(compo,CP_TEMPS,ti, tf, npas, 'rext', 6,
141 IF( info.NE. CPOK )GO TO 9000
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
147 c calcul du second membre
155 b(npt)=b(npt)+text(j)/(r/2+rflu(j))
159 b(npt)=b(npt)+puiss(npt)
162 c resolution du systeme lineaire
164 call dPBTRs( 'L' , N,kd, nrhs,A,LDAB,b,ldb,INFO )
175 write(6,*)"SOLIDE: temps=",tf
178 c ecriture de la temperature a t=tf
180 CALL cpeRE(compo,CP_TEMPS, tf, npas, 'temperature', 24,
182 IF( info.NE. CPOK )GO TO 9000
184 c ecriture de la temperature de paroi a t=tf
186 CALL cpeRE(compo,CP_TEMPS, tf, npas, 'tparoi', 6,
188 IF( info.NE. CPOK )GO TO 9000
190 c ecriture de la resistance de bord a t=tf
192 CALL cpeRE(compo,CP_TEMPS, tf, npas, 'rparoi', 6,
194 IF( info.NE. CPOK )GO TO 9000
200 CALL cpfin(compo,CP_ARRET, info)
203 SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
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
210 * .. Scalar Arguments ..
212 INTEGER INFO, KD, LDAB, N
214 * .. Array Arguments ..
215 DOUBLE PRECISION AB( LDAB, * )
221 * DPBTRF computes the Cholesky factorization of a real symmetric
222 * positive definite band matrix A.
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.
232 * UPLO (input) CHARACTER*1
233 * = 'U': Upper triangle of A is stored;
234 * = 'L': Lower triangle of A is stored.
237 * The order of the matrix A. N >= 0.
240 * The number of superdiagonals of the matrix A if UPLO = 'U',
241 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
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
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).
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.
255 * LDAB (input) INTEGER
256 * The leading dimension of the array AB. LDAB >= KD+1.
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
268 * The band storage scheme is illustrated by the following example, when
269 * N = 6, KD = 2, and UPLO = 'U':
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
277 * Similarly, if UPLO = 'L' the format of A is as follows:
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 * *
285 * Array elements marked * are not used by the routine.
288 * Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
290 * =====================================================================
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 )
298 * .. Local Scalars ..
299 INTEGER I, I2, I3, IB, II, J, JJ, NB
302 DOUBLE PRECISION WORK( LDWORK, NBMAX )
304 * .. External Functions ..
307 EXTERNAL LSAME, ILAENV
309 * .. External Subroutines ..
310 EXTERNAL DGEMM, DPBTF2, DPOTF2, DSYRK, DTRSM, XERBLA
312 * .. Intrinsic Functions ..
315 * .. Executable Statements ..
317 * Test the input parameters.
320 IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND.
321 $ ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
323 ELSE IF( N.LT.0 ) THEN
325 ELSE IF( KD.LT.0 ) THEN
327 ELSE IF( LDAB.LT.KD+1 ) THEN
331 CALL XERBLA( 'DPBTRF', -INFO )
335 * Quick return if possible
340 * Determine the block size for this environment
342 NB = ILAENV( 1, 'DPBTRF', UPLO, N, KD, -1, -1 )
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.
347 NB = MIN( NB, NBMAX )
349 IF( NB.LE.1 .OR. NB.GT.KD ) THEN
353 CALL DPBTF2( UPLO, N, KD, AB, LDAB, INFO )
358 IF( LSAME( UPLO, 'U' ) ) THEN
360 * Compute the Cholesky factorization of a symmetric band
361 * matrix, given the upper triangle of the matrix in band
364 * Zero the upper triangle of the work array.
372 * Process the band matrix one diagonal block at a time.
375 IB = MIN( NB, N-I+1 )
377 * Factorize the diagonal block
379 CALL DPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II )
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:
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.
400 I2 = MIN( KD-IB, N-I-IB+1 )
401 I3 = MIN( IB, N-I-KD+1 )
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 )
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 )
420 * Copy the lower triangle of A13 into the work array.
424 WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 )
428 * Update A13 (in the work array).
430 CALL DTRSM( 'Left', 'Upper', 'Transpose',
431 $ 'Non-unit', IB, I3, ONE, AB( KD+1, I ),
432 $ LDAB-1, WORK, LDWORK )
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 )
444 CALL DSYRK( 'Upper', 'Transpose', I3, IB, -ONE,
445 $ WORK, LDWORK, ONE, AB( KD+1, I+KD ),
448 * Copy the lower triangle of A13 back into place.
452 AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ )
460 * Compute the Cholesky factorization of a symmetric band
461 * matrix, given the lower triangle of the matrix in band
464 * Zero the lower triangle of the work array.
472 * Process the band matrix one diagonal block at a time.
475 IB = MIN( NB, N-I+1 )
477 * Factorize the diagonal block
479 CALL DPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II )
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:
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.
500 I2 = MIN( KD-IB, N-I-IB+1 )
501 I3 = MIN( IB, N-I-KD+1 )
507 CALL DTRSM( 'Right', 'Lower', 'Transpose',
508 $ 'Non-unit', I2, IB, ONE, AB( 1, I ),
509 $ LDAB-1, AB( 1+IB, I ), LDAB-1 )
513 CALL DSYRK( 'Lower', 'No Transpose', I2, IB, -ONE,
514 $ AB( 1+IB, I ), LDAB-1, ONE,
515 $ AB( 1, I+IB ), LDAB-1 )
520 * Copy the upper triangle of A31 into the work array.
523 DO 100 II = 1, MIN( JJ, I3 )
524 WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 )
528 * Update A31 (in the work array).
530 CALL DTRSM( 'Right', 'Lower', 'Transpose',
531 $ 'Non-unit', I3, IB, ONE, AB( 1, I ),
532 $ LDAB-1, WORK, LDWORK )
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 )
544 CALL DSYRK( 'Lower', 'No Transpose', I3, IB, -ONE,
545 $ WORK, LDWORK, ONE, AB( 1, I+KD ),
548 * Copy the upper triangle of A31 back into place.
551 DO 120 II = 1, MIN( JJ, I3 )
552 AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ )
568 SUBROUTINE DPBTRS( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO )
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
575 * .. Scalar Arguments ..
577 INTEGER INFO, KD, LDAB, LDB, N, NRHS
579 * .. Array Arguments ..
580 DOUBLE PRECISION AB( LDAB, * ), B( LDB, * )
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.
593 * UPLO (input) CHARACTER*1
594 * = 'U': Upper triangular factor stored in AB;
595 * = 'L': Lower triangular factor stored in AB.
598 * The order of the matrix A. N >= 0.
601 * The number of superdiagonals of the matrix A if UPLO = 'U',
602 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
604 * NRHS (input) INTEGER
605 * The number of right hand sides, i.e., the number of columns
606 * of the matrix B. NRHS >= 0.
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).
616 * LDAB (input) INTEGER
617 * The leading dimension of the array AB. LDAB >= KD+1.
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.
623 * LDB (input) INTEGER
624 * The leading dimension of the array B. LDB >= max(1,N).
626 * INFO (output) INTEGER
627 * = 0: successful exit
628 * < 0: if INFO = -i, the i-th argument had an illegal value
630 * =====================================================================
632 * .. Local Scalars ..
636 * .. External Functions ..
640 * .. External Subroutines ..
641 EXTERNAL DTBSV, XERBLA
643 * .. Intrinsic Functions ..
646 * .. Executable Statements ..
648 * Test the input parameters.
651 UPPER = LSAME( UPLO, 'U' )
652 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
654 ELSE IF( N.LT.0 ) THEN
656 ELSE IF( KD.LT.0 ) THEN
658 ELSE IF( NRHS.LT.0 ) THEN
660 ELSE IF( LDAB.LT.KD+1 ) THEN
662 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
666 CALL XERBLA( 'DPBTRS', -INFO )
670 * Quick return if possible
672 IF( N.EQ.0 .OR. NRHS.EQ.0 )
677 * Solve A*X = B where A = U'*U.
681 * Solve U'*X = B, overwriting B with X.
683 CALL DTBSV( 'Upper', 'Transpose', 'Non-unit', N, KD, AB,
684 $ LDAB, B( 1, J ), 1 )
686 * Solve U*X = B, overwriting B with X.
688 CALL DTBSV( 'Upper', 'No transpose', 'Non-unit', N, KD, AB,
689 $ LDAB, B( 1, J ), 1 )
693 * Solve A*X = B where A = L*L'.
697 * Solve L*X = B, overwriting B with X.
699 CALL DTBSV( 'Lower', 'No transpose', 'Non-unit', N, KD, AB,
700 $ LDAB, B( 1, J ), 1 )
702 * Solve L'*X = B, overwriting B with X.
704 CALL DTBSV( 'Lower', 'Transpose', 'Non-unit', N, KD, AB,
705 $ LDAB, B( 1, J ), 1 )
714 INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,
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
722 * .. Scalar Arguments ..
723 CHARACTER*( * ) NAME, OPTS
724 INTEGER ISPEC, N1, N2, N3, N4
730 * ILAENV is called from the LAPACK routines to choose problem-dependent
731 * parameters for the local environment. See ISPEC for a description of
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.
740 * This routine will not function correctly if it is converted to all
741 * lower case. Converting it to all upper case is allowed.
746 * ISPEC (input) INTEGER
747 * Specifies the parameter to be returned as the value of
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.
769 * NAME (input) CHARACTER*(*)
770 * The name of the calling subroutine, in either upper case or
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'.
783 * Problem dimensions for the subroutine NAME; these may not all
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.
793 * The following conventions have been used when calling ILAENV from the
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:
807 * NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
808 * IF( NB.LE.1 ) NB = MAX( 1, N )
810 * =====================================================================
812 * .. Local Scalars ..
818 INTEGER I, IC, IZ, NB, NBMIN, NX
820 * .. Intrinsic Functions ..
821 INTRINSIC CHAR, ICHAR, INT, MIN, REAL
823 * .. Executable Statements ..
825 GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC
827 * Invalid value for ISPEC
834 * Convert NAME to upper case if the first character is lower case.
838 IC = ICHAR( SUBNAM( 1:1 ) )
840 IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN
842 * ASCII character set
844 IF( IC.GE.97 .AND. IC.LE.122 ) THEN
845 SUBNAM( 1:1 ) = CHAR( IC-32 )
847 IC = ICHAR( SUBNAM( I:I ) )
848 IF( IC.GE.97 .AND. IC.LE.122 )
849 $ SUBNAM( I:I ) = CHAR( IC-32 )
853 ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN
855 * EBCDIC character set
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 )
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 )
870 ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN
872 * Prime machines: ASCII+128
874 IF( IC.GE.225 .AND. IC.LE.250 ) THEN
875 SUBNAM( 1:1 ) = CHAR( IC-32 )
877 IC = ICHAR( SUBNAM( I:I ) )
878 IF( IC.GE.225 .AND. IC.LE.250 )
879 $ SUBNAM( I:I ) = CHAR( IC-32 )
885 SNAME = C1.EQ.'S' .OR. C1.EQ.'D'
886 CNAME = C1.EQ.'C' .OR. C1.EQ.'Z'
887 IF( .NOT.( CNAME .OR. SNAME ) )
893 GO TO ( 110, 200, 300 ) ISPEC
897 * ISPEC = 1: block size
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.
905 IF( C2.EQ.'GE' ) THEN
906 IF( C3.EQ.'TRF' ) THEN
912 ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
919 ELSE IF( C3.EQ.'HRD' ) THEN
925 ELSE IF( C3.EQ.'BRD' ) THEN
931 ELSE IF( C3.EQ.'TRI' ) THEN
938 ELSE IF( C2.EQ.'PO' ) THEN
939 IF( C3.EQ.'TRF' ) THEN
946 ELSE IF( C2.EQ.'SY' ) THEN
947 IF( C3.EQ.'TRF' ) THEN
953 ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
955 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN
958 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
959 IF( C3.EQ.'TRF' ) THEN
961 ELSE IF( C3.EQ.'TRD' ) THEN
963 ELSE IF( C3.EQ.'GST' ) THEN
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.
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.
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.
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.
994 ELSE IF( C2.EQ.'GB' ) THEN
995 IF( C3.EQ.'TRF' ) THEN
1010 ELSE IF( C2.EQ.'PB' ) THEN
1011 IF( C3.EQ.'TRF' ) THEN
1026 ELSE IF( C2.EQ.'TR' ) THEN
1027 IF( C3.EQ.'TRI' ) THEN
1034 ELSE IF( C2.EQ.'LA' ) THEN
1035 IF( C3.EQ.'UUM' ) THEN
1042 ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
1043 IF( C3.EQ.'EBZ' ) THEN
1052 * ISPEC = 2: minimum block size
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
1063 ELSE IF( C3.EQ.'HRD' ) THEN
1069 ELSE IF( C3.EQ.'BRD' ) THEN
1075 ELSE IF( C3.EQ.'TRI' ) THEN
1082 ELSE IF( C2.EQ.'SY' ) THEN
1083 IF( C3.EQ.'TRF' ) THEN
1089 ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
1092 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
1093 IF( C3.EQ.'TRD' ) THEN
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.
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.
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.
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.
1130 * ISPEC = 3: crossover point
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
1141 ELSE IF( C3.EQ.'HRD' ) THEN
1147 ELSE IF( C3.EQ.'BRD' ) THEN
1154 ELSE IF( C2.EQ.'SY' ) THEN
1155 IF( SNAME .AND. C3.EQ.'TRD' ) THEN
1158 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
1159 IF( C3.EQ.'TRD' ) THEN
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.
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.
1184 * ISPEC = 4: number of shifts (used by xHSEQR)
1191 * ISPEC = 5: minimum column dimension (not used)
1198 * ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD)
1200 ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 )
1205 * ISPEC = 7: number of processors (not used)
1212 * ISPEC = 8: crossover point for multishift (used by xHSEQR)
1220 SUBROUTINE DPBTF2( UPLO, N, KD, AB, LDAB, INFO )
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
1227 * .. Scalar Arguments ..
1229 INTEGER INFO, KD, LDAB, N
1231 * .. Array Arguments ..
1232 DOUBLE PRECISION AB( LDAB, * )
1238 * DPBTF2 computes the Cholesky factorization of a real symmetric
1239 * positive definite band matrix A.
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.
1247 * This is the unblocked version of the algorithm, calling Level 2 BLAS.
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
1259 * The order of the matrix A. N >= 0.
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.
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
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).
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.
1277 * LDAB (input) INTEGER
1278 * The leading dimension of the array AB. LDAB >= KD+1.
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
1290 * The band storage scheme is illustrated by the following example, when
1291 * N = 6, KD = 2, and UPLO = 'U':
1293 * On entry: On exit:
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
1299 * Similarly, if UPLO = 'L' the format of A is as follows:
1301 * On entry: On exit:
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 * *
1307 * Array elements marked * are not used by the routine.
1309 * =====================================================================
1312 DOUBLE PRECISION ONE, ZERO
1313 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
1315 * .. Local Scalars ..
1318 DOUBLE PRECISION AJJ
1320 * .. External Functions ..
1324 * .. External Subroutines ..
1325 EXTERNAL DSCAL, DSYR, XERBLA
1327 * .. Intrinsic Functions ..
1328 INTRINSIC MAX, MIN, SQRT
1330 * .. Executable Statements ..
1332 * Test the input parameters.
1335 UPPER = LSAME( UPLO, 'U' )
1336 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
1338 ELSE IF( N.LT.0 ) THEN
1340 ELSE IF( KD.LT.0 ) THEN
1342 ELSE IF( LDAB.LT.KD+1 ) THEN
1345 IF( INFO.NE.0 ) THEN
1346 CALL XERBLA( 'DPBTF2', -INFO )
1350 * Quick return if possible
1355 KLD = MAX( 1, LDAB-1 )
1359 * Compute the Cholesky factorization A = U'*U.
1363 * Compute U(J,J) and test for non-positive-definiteness.
1371 * Compute elements J+1:J+KN of row J and update the
1372 * trailing submatrix within the band.
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 )
1383 * Compute the Cholesky factorization A = L*L'.
1387 * Compute L(J,J) and test for non-positive-definiteness.
1395 * Compute elements J+1:J+KN of column J and update the
1396 * trailing submatrix within the band.
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 )
1415 LOGICAL FUNCTION LSAME( CA, CB )
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
1422 * .. Scalar Arguments ..
1429 * LSAME returns .TRUE. if CA is the same letter as CB regardless of
1435 * CA (input) CHARACTER*1
1436 * CB (input) CHARACTER*1
1437 * CA and CB specify the single characters to be compared.
1439 * =====================================================================
1441 * .. Intrinsic Functions ..
1444 * .. Local Scalars ..
1445 INTEGER INTA, INTB, ZCODE
1447 * .. Executable Statements ..
1449 * Test if the characters are equal
1455 * Now test for equivalence if both characters are alphabetic.
1457 ZCODE = ICHAR( 'Z' )
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.
1467 IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN
1469 * ASCII is assumed - ZCODE is the ASCII code of either lower or
1472 IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32
1473 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32
1475 ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN
1477 * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
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
1487 ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN
1489 * ASCII is assumed, on Prime machines - ZCODE is the ASCII code
1490 * plus 128 of either lower or upper case 'Z'.
1492 IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32
1493 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32
1495 LSAME = INTA.EQ.INTB
1502 SUBROUTINE XERBLA( SRNAME, INFO )
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
1509 * .. Scalar Arguments ..
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.
1521 * Installers may consider modifying the STOP statement in order to
1522 * call system-specific exception-handling facilities.
1527 * SRNAME (input) CHARACTER*6
1528 * The name of the routine which called XERBLA.
1530 * INFO (input) INTEGER
1531 * The position of the invalid parameter in the parameter list
1532 * of the calling routine.
1534 * =====================================================================
1536 * .. Executable Statements ..
1538 WRITE( *, FMT = 9999 )SRNAME, INFO
1542 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ',
1543 $ 'an illegal value' )
1548 SUBROUTINE DPOTF2( UPLO, N, A, LDA, INFO )
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
1555 * .. Scalar Arguments ..
1557 INTEGER INFO, LDA, N
1559 * .. Array Arguments ..
1560 DOUBLE PRECISION A( LDA, * )
1566 * DPOTF2 computes the Cholesky factorization of a real symmetric
1567 * positive definite matrix A.
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.
1574 * This is the unblocked version of the algorithm, calling Level 2 BLAS.
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
1586 * The order of the matrix A. N >= 0.
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.
1597 * On exit, if INFO = 0, the factor U or L from the Cholesky
1598 * factorization A = U'*U or A = L*L'.
1600 * LDA (input) INTEGER
1601 * The leading dimension of the array A. LDA >= max(1,N).
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
1610 * =====================================================================
1613 DOUBLE PRECISION ONE, ZERO
1614 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
1616 * .. Local Scalars ..
1619 DOUBLE PRECISION AJJ
1621 * .. External Functions ..
1623 DOUBLE PRECISION DDOT
1624 EXTERNAL LSAME, DDOT
1626 * .. External Subroutines ..
1627 EXTERNAL DGEMV, DSCAL, XERBLA
1629 * .. Intrinsic Functions ..
1632 * .. Executable Statements ..
1634 * Test the input parameters.
1637 UPPER = LSAME( UPLO, 'U' )
1638 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
1640 ELSE IF( N.LT.0 ) THEN
1642 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
1645 IF( INFO.NE.0 ) THEN
1646 CALL XERBLA( 'DPOTF2', -INFO )
1650 * Quick return if possible
1657 * Compute the Cholesky factorization A = U'*U.
1661 * Compute U(J,J) and test for non-positive-definiteness.
1663 AJJ = A( J, J ) - DDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 )
1664 IF( AJJ.LE.ZERO ) THEN
1671 * Compute elements J+1:N of row J.
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 )
1681 * Compute the Cholesky factorization A = L*L'.
1685 * Compute L(J,J) and test for non-positive-definiteness.
1687 AJJ = A( J, J ) - DDOT( J-1, A( J, 1 ), LDA, A( J, 1 ),
1689 IF( AJJ.LE.ZERO ) THEN
1696 * Compute elements J+1:N of column J.
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 )
1716 subroutine dscal(n,da,dx,incx)
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(*)
1724 double precision da,dx(*)
1725 integer i,incx,m,mp1,n,nincx
1727 if( n.le.0 .or. incx.le.0 )return
1728 if(incx.eq.1)go to 20
1730 c code for increment not equal to 1
1733 do 10 i = 1,nincx,incx
1738 c code for increment equal to 1
1744 if( m .eq. 0 ) go to 40
1748 if( n .lt. 5 ) return
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)
1759 double precision function ddot(n,dx,incx,dy,incy)
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(*)
1766 double precision dx(*),dy(*),dtemp
1767 integer i,incx,incy,ix,iy,m,mp1,n
1772 if(incx.eq.1.and.incy.eq.1)go to 20
1774 c code for unequal increments or equal increments
1779 if(incx.lt.0)ix = (-n+1)*incx + 1
1780 if(incy.lt.0)iy = (-n+1)*incy + 1
1782 dtemp = dtemp + dx(ix)*dy(iy)
1789 c code for both increments equal to 1
1795 if( m .eq. 0 ) go to 40
1797 dtemp = dtemp + dx(i)*dy(i)
1799 if( n .lt. 5 ) go to 60
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)
1808 SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
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, * )
1821 * DTRSM solves one of the matrix equations
1823 * op( A )*X = alpha*B, or X*op( A ) = alpha*B,
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
1828 * op( A ) = A or op( A ) = A'.
1830 * The matrix X is overwritten on B.
1835 * SIDE - CHARACTER*1.
1836 * On entry, SIDE specifies whether op( A ) appears on the left
1837 * or right of X as follows:
1839 * SIDE = 'L' or 'l' op( A )*X = alpha*B.
1841 * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
1843 * Unchanged on exit.
1845 * UPLO - CHARACTER*1.
1846 * On entry, UPLO specifies whether the matrix A is an upper or
1847 * lower triangular matrix as follows:
1849 * UPLO = 'U' or 'u' A is an upper triangular matrix.
1851 * UPLO = 'L' or 'l' A is a lower triangular matrix.
1853 * Unchanged on exit.
1855 * TRANSA - CHARACTER*1.
1856 * On entry, TRANSA specifies the form of op( A ) to be used in
1857 * the matrix multiplication as follows:
1859 * TRANSA = 'N' or 'n' op( A ) = A.
1861 * TRANSA = 'T' or 't' op( A ) = A'.
1863 * TRANSA = 'C' or 'c' op( A ) = A'.
1865 * Unchanged on exit.
1867 * DIAG - CHARACTER*1.
1868 * On entry, DIAG specifies whether or not A is unit triangular
1871 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
1873 * DIAG = 'N' or 'n' A is not assumed to be unit
1876 * Unchanged on exit.
1879 * On entry, M specifies the number of rows of B. M must be at
1881 * Unchanged on exit.
1884 * On entry, N specifies the number of columns of B. N must be
1886 * Unchanged on exit.
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
1892 * Unchanged on exit.
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.
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.
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.
1921 * On entry, LDB specifies the first dimension of B as declared
1922 * in the calling (sub) program. LDB must be at least
1924 * Unchanged on exit.
1927 * Level 3 Blas routine.
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.
1937 * .. External Functions ..
1940 * .. External Subroutines ..
1942 * .. Intrinsic Functions ..
1944 * .. Local Scalars ..
1945 LOGICAL LSIDE, NOUNIT, UPPER
1946 INTEGER I, INFO, J, K, NROWA
1947 DOUBLE PRECISION TEMP
1949 DOUBLE PRECISION ONE , ZERO
1950 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
1952 * .. Executable Statements ..
1954 * Test the input parameters.
1956 LSIDE = LSAME( SIDE , 'L' )
1962 NOUNIT = LSAME( DIAG , 'N' )
1963 UPPER = LSAME( UPLO , 'U' )
1966 IF( ( .NOT.LSIDE ).AND.
1967 $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
1969 ELSE IF( ( .NOT.UPPER ).AND.
1970 $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
1972 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
1973 $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
1974 $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
1976 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
1977 $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
1979 ELSE IF( M .LT.0 )THEN
1981 ELSE IF( N .LT.0 )THEN
1983 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
1985 ELSE IF( LDB.LT.MAX( 1, M ) )THEN
1989 CALL XERBLA( 'DTRSM ', INFO )
1993 * Quick return if possible.
1998 * And when alpha.eq.zero.
2000 IF( ALPHA.EQ.ZERO )THEN
2009 * Start the operations.
2012 IF( LSAME( TRANSA, 'N' ) )THEN
2014 * Form B := alpha*inv( A )*B.
2018 IF( ALPHA.NE.ONE )THEN
2020 B( I, J ) = ALPHA*B( I, J )
2024 IF( B( K, J ).NE.ZERO )THEN
2026 $ B( K, J ) = B( K, J )/A( K, K )
2028 B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
2035 IF( ALPHA.NE.ONE )THEN
2037 B( I, J ) = ALPHA*B( I, J )
2041 IF( B( K, J ).NE.ZERO )THEN
2043 $ B( K, J ) = B( K, J )/A( K, K )
2045 B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
2053 * Form B := alpha*inv( A' )*B.
2058 TEMP = ALPHA*B( I, J )
2059 DO 110, K = 1, I - 1
2060 TEMP = TEMP - A( K, I )*B( K, J )
2063 $ TEMP = TEMP/A( I, I )
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 )
2075 $ TEMP = TEMP/A( I, I )
2082 IF( LSAME( TRANSA, 'N' ) )THEN
2084 * Form B := alpha*B*inv( A ).
2088 IF( ALPHA.NE.ONE )THEN
2090 B( I, J ) = ALPHA*B( I, J )
2093 DO 190, K = 1, J - 1
2094 IF( A( K, J ).NE.ZERO )THEN
2096 B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
2101 TEMP = ONE/A( J, J )
2103 B( I, J ) = TEMP*B( I, J )
2108 DO 260, J = N, 1, -1
2109 IF( ALPHA.NE.ONE )THEN
2111 B( I, J ) = ALPHA*B( I, J )
2114 DO 240, K = J + 1, N
2115 IF( A( K, J ).NE.ZERO )THEN
2117 B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
2122 TEMP = ONE/A( J, J )
2124 B( I, J ) = TEMP*B( I, J )
2131 * Form B := alpha*B*inv( A' ).
2134 DO 310, K = N, 1, -1
2136 TEMP = ONE/A( K, K )
2138 B( I, K ) = TEMP*B( I, K )
2141 DO 290, J = 1, K - 1
2142 IF( A( J, K ).NE.ZERO )THEN
2145 B( I, J ) = B( I, J ) - TEMP*B( I, K )
2149 IF( ALPHA.NE.ONE )THEN
2151 B( I, K ) = ALPHA*B( I, K )
2158 TEMP = ONE/A( K, K )
2160 B( I, K ) = TEMP*B( I, K )
2163 DO 340, J = K + 1, N
2164 IF( A( J, K ).NE.ZERO )THEN
2167 B( I, J ) = B( I, J ) - TEMP*B( I, K )
2171 IF( ALPHA.NE.ONE )THEN
2173 B( I, K ) = ALPHA*B( I, K )
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( * )
2197 * DTBSV solves one of the systems of equations
2199 * A*x = b, or A'*x = b,
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 )
2205 * No test for singularity or near-singularity is included in this
2206 * routine. Such tests must be performed before calling this routine.
2211 * UPLO - CHARACTER*1.
2212 * On entry, UPLO specifies whether the matrix is an upper or
2213 * lower triangular matrix as follows:
2215 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2217 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2219 * Unchanged on exit.
2221 * TRANS - CHARACTER*1.
2222 * On entry, TRANS specifies the equations to be solved as
2225 * TRANS = 'N' or 'n' A*x = b.
2227 * TRANS = 'T' or 't' A'*x = b.
2229 * TRANS = 'C' or 'c' A'*x = b.
2231 * Unchanged on exit.
2233 * DIAG - CHARACTER*1.
2234 * On entry, DIAG specifies whether or not A is unit
2235 * triangular as follows:
2237 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2239 * DIAG = 'N' or 'n' A is not assumed to be unit
2242 * Unchanged on exit.
2245 * On entry, N specifies the order of the matrix A.
2246 * N must be at least zero.
2247 * Unchanged on exit.
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.
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
2271 * DO 10, I = MAX( 1, J - K ), J
2272 * A( M + I, J ) = matrix( I, J )
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
2289 * DO 10, I = J, MIN( N, J + K )
2290 * A( M + I, J ) = matrix( I, J )
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.
2300 * On entry, LDA specifies the first dimension of A as declared
2301 * in the calling (sub) program. LDA must be at least
2303 * Unchanged on exit.
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.
2312 * On entry, INCX specifies the increment for the elements of
2313 * X. INCX must not be zero.
2314 * Unchanged on exit.
2317 * Level 2 Blas routine.
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.
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
2333 * .. External Functions ..
2336 * .. External Subroutines ..
2338 * .. Intrinsic Functions ..
2341 * .. Executable Statements ..
2343 * Test the input parameters.
2346 IF ( .NOT.LSAME( UPLO , 'U' ).AND.
2347 $ .NOT.LSAME( UPLO , 'L' ) )THEN
2349 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
2350 $ .NOT.LSAME( TRANS, 'T' ).AND.
2351 $ .NOT.LSAME( TRANS, 'C' ) )THEN
2353 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
2354 $ .NOT.LSAME( DIAG , 'N' ) )THEN
2356 ELSE IF( N.LT.0 )THEN
2358 ELSE IF( K.LT.0 )THEN
2360 ELSE IF( LDA.LT.( K + 1 ) )THEN
2362 ELSE IF( INCX.EQ.0 )THEN
2366 CALL XERBLA( 'DTBSV ', INFO )
2370 * Quick return if possible.
2375 NOUNIT = LSAME( DIAG, 'N' )
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.
2381 KX = 1 - ( N - 1 )*INCX
2382 ELSE IF( INCX.NE.1 )THEN
2386 * Start the operations. In this version the elements of A are
2387 * accessed by sequentially with one pass through A.
2389 IF( LSAME( TRANS, 'N' ) )THEN
2391 * Form x := inv( A )*x.
2393 IF( LSAME( UPLO, 'U' ) )THEN
2397 IF( X( J ).NE.ZERO )THEN
2400 $ X( J ) = X( J )/A( KPLUS1, J )
2402 DO 10, I = J - 1, MAX( 1, J - K ), -1
2403 X( I ) = X( I ) - TEMP*A( L + I, J )
2408 KX = KX + ( N - 1 )*INCX
2412 IF( X( JX ).NE.ZERO )THEN
2416 $ X( JX ) = X( JX )/A( KPLUS1, J )
2418 DO 30, I = J - 1, MAX( 1, J - K ), -1
2419 X( IX ) = X( IX ) - TEMP*A( L + I, J )
2429 IF( X( J ).NE.ZERO )THEN
2432 $ X( J ) = X( J )/A( 1, J )
2434 DO 50, I = J + 1, MIN( N, J + K )
2435 X( I ) = X( I ) - TEMP*A( L + I, J )
2443 IF( X( JX ).NE.ZERO )THEN
2447 $ X( JX ) = X( JX )/A( 1, J )
2449 DO 70, I = J + 1, MIN( N, J + K )
2450 X( IX ) = X( IX ) - TEMP*A( L + I, J )
2460 * Form x := inv( A')*x.
2462 IF( LSAME( UPLO, 'U' ) )THEN
2468 DO 90, I = MAX( 1, J - K ), J - 1
2469 TEMP = TEMP - A( L + I, J )*X( I )
2472 $ TEMP = TEMP/A( KPLUS1, J )
2481 DO 110, I = MAX( 1, J - K ), J - 1
2482 TEMP = TEMP - A( L + I, J )*X( IX )
2486 $ TEMP = TEMP/A( KPLUS1, J )
2495 DO 140, J = N, 1, -1
2498 DO 130, I = MIN( N, J + K ), J + 1, -1
2499 TEMP = TEMP - A( L + I, J )*X( I )
2502 $ TEMP = TEMP/A( 1, J )
2506 KX = KX + ( N - 1 )*INCX
2508 DO 160, J = N, 1, -1
2512 DO 150, I = MIN( N, J + K ), J + 1, -1
2513 TEMP = TEMP - A( L + I, J )*X( IX )
2517 $ TEMP = TEMP/A( 1, J )
2520 IF( ( N - J ).GE.K )
2532 SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA )
2533 * .. Scalar Arguments ..
2534 DOUBLE PRECISION ALPHA
2535 INTEGER INCX, LDA, N
2537 * .. Array Arguments ..
2538 DOUBLE PRECISION A( LDA, * ), X( * )
2544 * DSYR performs the symmetric rank 1 operation
2546 * A := alpha*x*x' + A,
2548 * where alpha is a real scalar, x is an n element vector and A is an
2549 * n by n symmetric matrix.
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
2559 * UPLO = 'U' or 'u' Only the upper triangular part of A
2560 * is to be referenced.
2562 * UPLO = 'L' or 'l' Only the lower triangular part of A
2563 * is to be referenced.
2565 * Unchanged on exit.
2568 * On entry, N specifies the order of the matrix A.
2569 * N must be at least zero.
2570 * Unchanged on exit.
2572 * ALPHA - DOUBLE PRECISION.
2573 * On entry, ALPHA specifies the scalar alpha.
2574 * Unchanged on exit.
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
2580 * Unchanged on exit.
2583 * On entry, INCX specifies the increment for the elements of
2584 * X. INCX must not be zero.
2585 * Unchanged on exit.
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.
2602 * On entry, LDA specifies the first dimension of A as declared
2603 * in the calling (sub) program. LDA must be at least
2605 * Unchanged on exit.
2608 * Level 2 Blas routine.
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.
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 ..
2626 * .. External Subroutines ..
2628 * .. Intrinsic Functions ..
2631 * .. Executable Statements ..
2633 * Test the input parameters.
2636 IF ( .NOT.LSAME( UPLO, 'U' ).AND.
2637 $ .NOT.LSAME( UPLO, 'L' ) )THEN
2639 ELSE IF( N.LT.0 )THEN
2641 ELSE IF( INCX.EQ.0 )THEN
2643 ELSE IF( LDA.LT.MAX( 1, N ) )THEN
2647 CALL XERBLA( 'DSYR ', INFO )
2651 * Quick return if possible.
2653 IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
2656 * Set the start point in X if the increment is not unity.
2659 KX = 1 - ( N - 1 )*INCX
2660 ELSE IF( INCX.NE.1 )THEN
2664 * Start the operations. In this version the elements of A are
2665 * accessed sequentially with one pass through the triangular part
2668 IF( LSAME( UPLO, 'U' ) )THEN
2670 * Form A when A is stored in upper triangle.
2674 IF( X( J ).NE.ZERO )THEN
2677 A( I, J ) = A( I, J ) + X( I )*TEMP
2684 IF( X( JX ).NE.ZERO )THEN
2685 TEMP = ALPHA*X( JX )
2688 A( I, J ) = A( I, J ) + X( IX )*TEMP
2697 * Form A when A is stored in lower triangle.
2701 IF( X( J ).NE.ZERO )THEN
2704 A( I, J ) = A( I, J ) + X( I )*TEMP
2711 IF( X( JX ).NE.ZERO )THEN
2712 TEMP = ALPHA*X( JX )
2715 A( I, J ) = A( I, J ) + X( IX )*TEMP
2729 SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
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, * )
2742 * DSYRK performs one of the symmetric rank k operations
2744 * C := alpha*A*A' + beta*C,
2748 * C := alpha*A'*A + beta*C,
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.
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
2762 * UPLO = 'U' or 'u' Only the upper triangular part of C
2763 * is to be referenced.
2765 * UPLO = 'L' or 'l' Only the lower triangular part of C
2766 * is to be referenced.
2768 * Unchanged on exit.
2770 * TRANS - CHARACTER*1.
2771 * On entry, TRANS specifies the operation to be performed as
2774 * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
2776 * TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
2778 * TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
2780 * Unchanged on exit.
2783 * On entry, N specifies the order of the matrix C. N must be
2785 * Unchanged on exit.
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.
2794 * ALPHA - DOUBLE PRECISION.
2795 * On entry, ALPHA specifies the scalar alpha.
2796 * Unchanged on exit.
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
2804 * Unchanged on exit.
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.
2813 * BETA - DOUBLE PRECISION.
2814 * On entry, BETA specifies the scalar beta.
2815 * Unchanged on exit.
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.
2832 * On entry, LDC specifies the first dimension of C as declared
2833 * in the calling (sub) program. LDC must be at least
2835 * Unchanged on exit.
2838 * Level 3 Blas routine.
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.
2847 * .. External Functions ..
2850 * .. External Subroutines ..
2852 * .. Intrinsic Functions ..
2854 * .. Local Scalars ..
2856 INTEGER I, INFO, J, L, NROWA
2857 DOUBLE PRECISION TEMP
2859 DOUBLE PRECISION ONE , ZERO
2860 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
2862 * .. Executable Statements ..
2864 * Test the input parameters.
2866 IF( LSAME( TRANS, 'N' ) )THEN
2871 UPPER = LSAME( UPLO, 'U' )
2874 IF( ( .NOT.UPPER ).AND.
2875 $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
2877 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
2878 $ ( .NOT.LSAME( TRANS, 'T' ) ).AND.
2879 $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN
2881 ELSE IF( N .LT.0 )THEN
2883 ELSE IF( K .LT.0 )THEN
2885 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
2887 ELSE IF( LDC.LT.MAX( 1, N ) )THEN
2891 CALL XERBLA( 'DSYRK ', INFO )
2895 * Quick return if possible.
2898 $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
2901 * And when alpha.eq.zero.
2903 IF( ALPHA.EQ.ZERO )THEN
2905 IF( BETA.EQ.ZERO )THEN
2914 C( I, J ) = BETA*C( I, J )
2919 IF( BETA.EQ.ZERO )THEN
2928 C( I, J ) = BETA*C( I, J )
2936 * Start the operations.
2938 IF( LSAME( TRANS, 'N' ) )THEN
2940 * Form C := alpha*A*A' + beta*C.
2944 IF( BETA.EQ.ZERO )THEN
2948 ELSE IF( BETA.NE.ONE )THEN
2950 C( I, J ) = BETA*C( I, J )
2954 IF( A( J, L ).NE.ZERO )THEN
2955 TEMP = ALPHA*A( J, L )
2957 C( I, J ) = C( I, J ) + TEMP*A( I, L )
2964 IF( BETA.EQ.ZERO )THEN
2968 ELSE IF( BETA.NE.ONE )THEN
2970 C( I, J ) = BETA*C( I, J )
2974 IF( A( J, L ).NE.ZERO )THEN
2975 TEMP = ALPHA*A( J, L )
2977 C( I, J ) = C( I, J ) + TEMP*A( I, L )
2985 * Form C := alpha*A'*A + beta*C.
2992 TEMP = TEMP + A( L, I )*A( L, J )
2994 IF( BETA.EQ.ZERO )THEN
2995 C( I, J ) = ALPHA*TEMP
2997 C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
3006 TEMP = TEMP + A( L, I )*A( L, J )
3008 IF( BETA.EQ.ZERO )THEN
3009 C( I, J ) = ALPHA*TEMP
3011 C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
3023 SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
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, * )
3036 * DGEMM performs one of the matrix-matrix operations
3038 * C := alpha*op( A )*op( B ) + beta*C,
3040 * where op( X ) is one of
3042 * op( X ) = X or op( X ) = X',
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.
3050 * TRANSA - CHARACTER*1.
3051 * On entry, TRANSA specifies the form of op( A ) to be used in
3052 * the matrix multiplication as follows:
3054 * TRANSA = 'N' or 'n', op( A ) = A.
3056 * TRANSA = 'T' or 't', op( A ) = A'.
3058 * TRANSA = 'C' or 'c', op( A ) = A'.
3060 * Unchanged on exit.
3062 * TRANSB - CHARACTER*1.
3063 * On entry, TRANSB specifies the form of op( B ) to be used in
3064 * the matrix multiplication as follows:
3066 * TRANSB = 'N' or 'n', op( B ) = B.
3068 * TRANSB = 'T' or 't', op( B ) = B'.
3070 * TRANSB = 'C' or 'c', op( B ) = B'.
3072 * Unchanged on exit.
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.
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
3083 * Unchanged on exit.
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
3089 * Unchanged on exit.
3091 * ALPHA - DOUBLE PRECISION.
3092 * On entry, ALPHA specifies the scalar alpha.
3093 * Unchanged on exit.
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
3101 * Unchanged on exit.
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.
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
3116 * Unchanged on exit.
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.
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.
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 ).
3138 * On entry, LDC specifies the first dimension of C as declared
3139 * in the calling (sub) program. LDC must be at least
3141 * Unchanged on exit.
3144 * Level 3 Blas routine.
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.
3153 * .. External Functions ..
3156 * .. External Subroutines ..
3158 * .. Intrinsic Functions ..
3160 * .. Local Scalars ..
3162 INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB
3163 DOUBLE PRECISION TEMP
3165 DOUBLE PRECISION ONE , ZERO
3166 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
3168 * .. Executable Statements ..
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.
3174 NOTA = LSAME( TRANSA, 'N' )
3175 NOTB = LSAME( TRANSB, 'N' )
3189 * Test the input parameters.
3192 IF( ( .NOT.NOTA ).AND.
3193 $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
3194 $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN
3196 ELSE IF( ( .NOT.NOTB ).AND.
3197 $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
3198 $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN
3200 ELSE IF( M .LT.0 )THEN
3202 ELSE IF( N .LT.0 )THEN
3204 ELSE IF( K .LT.0 )THEN
3206 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
3208 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
3210 ELSE IF( LDC.LT.MAX( 1, M ) )THEN
3214 CALL XERBLA( 'DGEMM ', INFO )
3218 * Quick return if possible.
3220 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
3221 $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
3224 * And if alpha.eq.zero.
3226 IF( ALPHA.EQ.ZERO )THEN
3227 IF( BETA.EQ.ZERO )THEN
3236 C( I, J ) = BETA*C( I, J )
3243 * Start the operations.
3248 * Form C := alpha*A*B + beta*C.
3251 IF( BETA.EQ.ZERO )THEN
3255 ELSE IF( BETA.NE.ONE )THEN
3257 C( I, J ) = BETA*C( I, J )
3261 IF( B( L, J ).NE.ZERO )THEN
3262 TEMP = ALPHA*B( L, J )
3264 C( I, J ) = C( I, J ) + TEMP*A( I, L )
3271 * Form C := alpha*A'*B + beta*C
3277 TEMP = TEMP + A( L, I )*B( L, J )
3279 IF( BETA.EQ.ZERO )THEN
3280 C( I, J ) = ALPHA*TEMP
3282 C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
3290 * Form C := alpha*A*B' + beta*C
3293 IF( BETA.EQ.ZERO )THEN
3297 ELSE IF( BETA.NE.ONE )THEN
3299 C( I, J ) = BETA*C( I, J )
3303 IF( B( J, L ).NE.ZERO )THEN
3304 TEMP = ALPHA*B( J, L )
3306 C( I, J ) = C( I, J ) + TEMP*A( I, L )
3313 * Form C := alpha*A'*B' + beta*C
3319 TEMP = TEMP + A( L, I )*B( J, L )
3321 IF( BETA.EQ.ZERO )THEN
3322 C( I, J ) = ALPHA*TEMP
3324 C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
3336 SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
3338 * .. Scalar Arguments ..
3339 DOUBLE PRECISION ALPHA, BETA
3340 INTEGER INCX, INCY, LDA, M, N
3342 * .. Array Arguments ..
3343 DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )
3349 * DGEMV performs one of the matrix-vector operations
3351 * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
3353 * where alpha and beta are scalars, x and y are vectors and A is an
3359 * TRANS - CHARACTER*1.
3360 * On entry, TRANS specifies the operation to be performed as
3363 * TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
3365 * TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
3367 * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
3369 * Unchanged on exit.
3372 * On entry, M specifies the number of rows of the matrix A.
3373 * M must be at least zero.
3374 * Unchanged on exit.
3377 * On entry, N specifies the number of columns of the matrix A.
3378 * N must be at least zero.
3379 * Unchanged on exit.
3381 * ALPHA - DOUBLE PRECISION.
3382 * On entry, ALPHA specifies the scalar alpha.
3383 * Unchanged on exit.
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.
3391 * On entry, LDA specifies the first dimension of A as declared
3392 * in the calling (sub) program. LDA must be at least
3394 * Unchanged on exit.
3396 * X - DOUBLE PRECISION array of DIMENSION at least
3397 * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
3399 * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
3400 * Before entry, the incremented array X must contain the
3402 * Unchanged on exit.
3405 * On entry, INCX specifies the increment for the elements of
3406 * X. INCX must not be zero.
3407 * Unchanged on exit.
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.
3414 * Y - DOUBLE PRECISION array of DIMENSION at least
3415 * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
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
3423 * On entry, INCY specifies the increment for the elements of
3424 * Y. INCY must not be zero.
3425 * Unchanged on exit.
3428 * Level 2 Blas routine.
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.
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 ..
3446 * .. External Subroutines ..
3448 * .. Intrinsic Functions ..
3451 * .. Executable Statements ..
3453 * Test the input parameters.
3456 IF ( .NOT.LSAME( TRANS, 'N' ).AND.
3457 $ .NOT.LSAME( TRANS, 'T' ).AND.
3458 $ .NOT.LSAME( TRANS, 'C' ) )THEN
3460 ELSE IF( M.LT.0 )THEN
3462 ELSE IF( N.LT.0 )THEN
3464 ELSE IF( LDA.LT.MAX( 1, M ) )THEN
3466 ELSE IF( INCX.EQ.0 )THEN
3468 ELSE IF( INCY.EQ.0 )THEN
3472 CALL XERBLA( 'DGEMV ', INFO )
3476 * Quick return if possible.
3478 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
3479 $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
3482 * Set LENX and LENY, the lengths of the vectors x and y, and set
3483 * up the start points in X and Y.
3485 IF( LSAME( TRANS, 'N' ) )THEN
3495 KX = 1 - ( LENX - 1 )*INCX
3500 KY = 1 - ( LENY - 1 )*INCY
3503 * Start the operations. In this version the elements of A are
3504 * accessed sequentially with one pass through A.
3506 * First form y := beta*y.
3508 IF( BETA.NE.ONE )THEN
3510 IF( BETA.EQ.ZERO )THEN
3516 Y( I ) = BETA*Y( I )
3521 IF( BETA.EQ.ZERO )THEN
3528 Y( IY ) = BETA*Y( IY )
3536 IF( LSAME( TRANS, 'N' ) )THEN
3538 * Form y := alpha*A*x + y.
3543 IF( X( JX ).NE.ZERO )THEN
3544 TEMP = ALPHA*X( JX )
3546 Y( I ) = Y( I ) + TEMP*A( I, J )
3553 IF( X( JX ).NE.ZERO )THEN
3554 TEMP = ALPHA*X( JX )
3557 Y( IY ) = Y( IY ) + TEMP*A( I, J )
3566 * Form y := alpha*A'*x + y.
3573 TEMP = TEMP + A( I, J )*X( I )
3575 Y( JY ) = Y( JY ) + ALPHA*TEMP
3583 TEMP = TEMP + A( I, J )*X( IX )
3586 Y( JY ) = Y( JY ) + ALPHA*TEMP
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)
3624 c construction de la matrice (laplacien)
3656 a(1,npt)=3./r+1./(r/2.+rext)
3669 a(1,i)=2./r+1./(r/2.+rext)
3673 a(1,i)=2./r+1./(r/2.+rext)
3678 c factorisation de la matrice
3683 call dPBTRF( 'L' , N, KD, A, LDAB, INFO )
3702 c initialisation de la temperature a iter=0
3704 CALL cpeRE(compo,CP_ITERATION, ti, iter, 'tempi', 24,
3706 IF( info.NE. CPOK )GO TO 9000
3708 c initialisation de la temperature de bord a iter=0.
3710 CALL cpeRE(compo,CP_ITERATION, ti, iter, 'tpi', 6,
3712 IF( info.NE. CPOK )GO TO 9000
3714 c boucle iterative infinie
3716 do while( iter .lt. 500)
3718 c lecture de la puissance combustible a iter
3720 CALL cplRE(compo,CP_ITERATION,ti, tf, iter, 'puissi', 24,
3721 & nval, puiss , info)
3722 IF( info.NE. CPOK )GO TO 9000
3724 c lecture de la temperature exterieure a iter
3726 CALL cplRE(compo,CP_ITERATION,ti, tf, iter, 'tfi', 6,
3727 & nval, text , info)
3728 IF( info.NE. CPOK )GO TO 9000
3730 c calcul du second membre
3738 b(npt)=b(npt)+text(j)/(r/2+rflu(j))
3742 c resolution du systeme lineaire
3744 call dPBTRs( 'L' , N,kd, nrhs,A,LDAB,b,ldb,INFO )
3756 c ecriture de la temperature a iter+1
3758 CALL cpeRE(compo,CP_ITERATION, ti, iter, 'tempi', 24,
3760 IF( info.NE. CPOK )GO TO 9000
3762 c ecriture de la temperature de paroi a iter+1
3764 CALL cpeRE(compo,CP_ITERATION, ti, iter, 'tpi', 6,
3766 IF( info.NE. CPOK )GO TO 9000
3768 c calcul de convergence
3772 conv1=conv1+(puiss(i)-puissn(i))**2
3773 conv1=conv1+(t(i)-tn(i))**2
3776 conv1=conv1+(tpar(i)-tparn(i))**2
3777 conv1=conv1+(text(i)-textn(i))**2
3781 if(conv1.lt.1.e-3) iconv=1
3783 c ecriture du flag de convergence iconv
3785 write(6,*)"SOLIDE:",iter,iconv
3787 CALL cpeEN(compo,CP_ITERATION,tf, iter , 'iconv', 1,
3789 IF( info.NE. CPOK )GO TO 9000
3791 if(iconv.eq.1)go to 9000
3805 CALL cpfin(compo,CP_ARRET, info)