1! 2! www.netlib.org/lapack/LICENSE.txt 3! 4! Copyright (c) 1992-2013 The University of Tennessee and The University 5! of Tennessee Research Foundation. All rights 6! reserved. 7! Copyright (c) 2000-2013 The University of California Berkeley. All 8! rights reserved. 9! Copyright (c) 2006-2013 The University of Colorado Denver. All rights 10! reserved. 11! 12! $COPYRIGHT$ 13! 14! Additional copyrights may follow 15! 16! $HEADER$ 17! 18! Redistribution and use in source and binary forms, with or without 19! modification, are permitted provided that the following conditions are 20! met: 21! 22! - Redistributions of source code must retain the above copyright 23! notice, this list of conditions and the following disclaimer. 24! 25! - Redistributions in binary form must reproduce the above copyright 26! notice, this list of conditions and the following disclaimer listed 27! in this license in the documentation and/or other materials 28! provided with the distribution. 29! 30! - Neither the name of the copyright holders nor the names of its 31! contributors may be used to endorse or promote products derived from 32! this software without specific prior written permission. 33! 34! The copyright holders provide no reassurances that the source code 35! provided does not infringe any patent, copyright, or any other 36! intellectual property rights of third parties. The copyright holders 37! disclaim any liability to any recipient for claims brought against 38! recipient by any third party for infringement of that parties 39! intellectual property rights. 40! 41! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 42! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 43! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 44! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 45! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 46! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 47! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 48! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 49! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 50! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 51! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 52! 53! subroutines to solve a set of linear equations with 54! a general real matrix 55! 56! 57 SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) 58* 59* -- LAPACK driver routine (version 3.0) -- 60* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 61* Courant Institute, Argonne National Lab, and Rice University 62* March 31, 1993 63* 64* .. Scalar Arguments .. 65 INTEGER INFO, LDA, LDB, N, NRHS 66* .. 67* .. Array Arguments .. 68 INTEGER IPIV( * ) 69 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 70* .. 71* 72* Purpose 73* ======= 74* 75* DGESV computes the solution to a real system of linear equations 76* A * X = B, 77* where A is an N-by-N matrix and X and B are N-by-NRHS matrices. 78* 79* The LU decomposition with partial pivoting and row interchanges is 80* used to factor A as 81* A = P * L * U, 82* where P is a permutation matrix, L is unit lower triangular, and U is 83* upper triangular. The factored form of A is then used to solve the 84* system of equations A * X = B. 85* 86* Arguments 87* ========= 88* 89* N (input) INTEGER 90* The number of linear equations, i.e., the order of the 91* matrix A. N >= 0. 92* 93* NRHS (input) INTEGER 94* The number of right hand sides, i.e., the number of columns 95* of the matrix B. NRHS >= 0. 96* 97* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 98* On entry, the N-by-N coefficient matrix A. 99* On exit, the factors L and U from the factorization 100* A = P*L*U; the unit diagonal elements of L are not stored. 101* 102* LDA (input) INTEGER 103* The leading dimension of the array A. LDA >= max(1,N). 104* 105* IPIV (output) INTEGER array, dimension (N) 106* The pivot indices that define the permutation matrix P; 107* row i of the matrix was interchanged with row IPIV(i). 108* 109* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 110* On entry, the N-by-NRHS matrix of right hand side matrix B. 111* On exit, if INFO = 0, the N-by-NRHS solution matrix X. 112* 113* LDB (input) INTEGER 114* The leading dimension of the array B. LDB >= max(1,N). 115* 116* INFO (output) INTEGER 117* = 0: successful exit 118* < 0: if INFO = -i, the i-th argument had an illegal value 119* > 0: if INFO = i, U(i,i) is exactly zero. The factorization 120* has been completed, but the factor U is exactly 121* singular, so the solution could not be computed. 122* 123* ===================================================================== 124* 125* .. External Subroutines .. 126 EXTERNAL DGETRF, DGETRS, XERBLA 127* .. 128* .. Intrinsic Functions .. 129 INTRINSIC MAX 130* .. 131* .. Executable Statements .. 132* 133* Test the input parameters. 134* 135 INFO = 0 136 IF( N.LT.0 ) THEN 137 INFO = -1 138 ELSE IF( NRHS.LT.0 ) THEN 139 INFO = -2 140 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 141 INFO = -4 142 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 143 INFO = -7 144 END IF 145 IF( INFO.NE.0 ) THEN 146 CALL XERBLA( 'DGESV ', -INFO ) 147 RETURN 148 END IF 149* 150* Compute the LU factorization of A. 151* 152 CALL DGETRF( N, N, A, LDA, IPIV, INFO ) 153 IF( INFO.EQ.0 ) THEN 154* 155* Solve the system A*X = B, overwriting B with X. 156* 157 CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, B, LDB, 158 $ INFO ) 159 END IF 160 RETURN 161* 162* End of DGESV 163* 164 END 165 SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) 166* 167* -- LAPACK routine (version 3.0) -- 168* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 169* Courant Institute, Argonne National Lab, and Rice University 170* June 30, 1992 171* 172* .. Scalar Arguments .. 173 INTEGER INFO, LDA, M, N 174* .. 175* .. Array Arguments .. 176 INTEGER IPIV( * ) 177 DOUBLE PRECISION A( LDA, * ) 178* .. 179* 180* Purpose 181* ======= 182* 183* DGETF2 computes an LU factorization of a general m-by-n matrix A 184* using partial pivoting with row interchanges. 185* 186* The factorization has the form 187* A = P * L * U 188* where P is a permutation matrix, L is lower triangular with unit 189* diagonal elements (lower trapezoidal if m > n), and U is upper 190* triangular (upper trapezoidal if m < n). 191* 192* This is the right-looking Level 2 BLAS version of the algorithm. 193* 194* Arguments 195* ========= 196* 197* M (input) INTEGER 198* The number of rows of the matrix A. M >= 0. 199* 200* N (input) INTEGER 201* The number of columns of the matrix A. N >= 0. 202* 203* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 204* On entry, the m by n matrix to be factored. 205* On exit, the factors L and U from the factorization 206* A = P*L*U; the unit diagonal elements of L are not stored. 207* 208* LDA (input) INTEGER 209* The leading dimension of the array A. LDA >= max(1,M). 210* 211* IPIV (output) INTEGER array, dimension (min(M,N)) 212* The pivot indices; for 1 <= i <= min(M,N), row i of the 213* matrix was interchanged with row IPIV(i). 214* 215* INFO (output) INTEGER 216* = 0: successful exit 217* < 0: if INFO = -k, the k-th argument had an illegal value 218* > 0: if INFO = k, U(k,k) is exactly zero. The factorization 219* has been completed, but the factor U is exactly 220* singular, and division by zero will occur if it is used 221* to solve a system of equations. 222* 223* ===================================================================== 224* 225* .. Parameters .. 226 DOUBLE PRECISION ONE, ZERO 227 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 228* .. 229* .. Local Scalars .. 230 INTEGER J, JP 231* .. 232* .. External Functions .. 233 INTEGER IDAMAX 234 EXTERNAL IDAMAX 235* .. 236* .. External Subroutines .. 237 EXTERNAL DGER, DSCAL, DSWAP, XERBLA 238* .. 239* .. Intrinsic Functions .. 240 INTRINSIC MAX, MIN 241* .. 242* .. Executable Statements .. 243* 244* Test the input parameters. 245* 246 INFO = 0 247 IF( M.LT.0 ) THEN 248 INFO = -1 249 ELSE IF( N.LT.0 ) THEN 250 INFO = -2 251 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 252 INFO = -4 253 END IF 254 IF( INFO.NE.0 ) THEN 255 CALL XERBLA( 'DGETF2', -INFO ) 256 RETURN 257 END IF 258* 259* Quick return if possible 260* 261 IF( M.EQ.0 .OR. N.EQ.0 ) 262 $ RETURN 263* 264 DO 10 J = 1, MIN( M, N ) 265* 266* Find pivot and test for singularity. 267* 268 JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 ) 269 IPIV( J ) = JP 270 IF( A( JP, J ).NE.ZERO ) THEN 271* 272* Apply the interchange to columns 1:N. 273* 274 IF( JP.NE.J ) 275 $ CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) 276* 277* Compute elements J+1:M of J-th column. 278* 279 IF( J.LT.M ) 280 $ CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) 281* 282 ELSE IF( INFO.EQ.0 ) THEN 283* 284 INFO = J 285 END IF 286* 287 IF( J.LT.MIN( M, N ) ) THEN 288* 289* Update trailing submatrix. 290* 291 CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA, 292 $ A( J+1, J+1 ), LDA ) 293 END IF 294 10 CONTINUE 295 RETURN 296* 297* End of DGETF2 298* 299 END 300 SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) 301* 302* -- LAPACK routine (version 3.0) -- 303* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 304* Courant Institute, Argonne National Lab, and Rice University 305* March 31, 1993 306* 307* .. Scalar Arguments .. 308 INTEGER INFO, LDA, M, N 309* .. 310* .. Array Arguments .. 311 INTEGER IPIV( * ) 312 DOUBLE PRECISION A( LDA, * ) 313* .. 314* 315* Purpose 316* ======= 317* 318* DGETRF computes an LU factorization of a general M-by-N matrix A 319* using partial pivoting with row interchanges. 320* 321* The factorization has the form 322* A = P * L * U 323* where P is a permutation matrix, L is lower triangular with unit 324* diagonal elements (lower trapezoidal if m > n), and U is upper 325* triangular (upper trapezoidal if m < n). 326* 327* This is the right-looking Level 3 BLAS version of the algorithm. 328* 329* Arguments 330* ========= 331* 332* M (input) INTEGER 333* The number of rows of the matrix A. M >= 0. 334* 335* N (input) INTEGER 336* The number of columns of the matrix A. N >= 0. 337* 338* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 339* On entry, the M-by-N matrix to be factored. 340* On exit, the factors L and U from the factorization 341* A = P*L*U; the unit diagonal elements of L are not stored. 342* 343* LDA (input) INTEGER 344* The leading dimension of the array A. LDA >= max(1,M). 345* 346* IPIV (output) INTEGER array, dimension (min(M,N)) 347* The pivot indices; for 1 <= i <= min(M,N), row i of the 348* matrix was interchanged with row IPIV(i). 349* 350* INFO (output) INTEGER 351* = 0: successful exit 352* < 0: if INFO = -i, the i-th argument had an illegal value 353* > 0: if INFO = i, U(i,i) is exactly zero. The factorization 354* has been completed, but the factor U is exactly 355* singular, and division by zero will occur if it is used 356* to solve a system of equations. 357* 358* ===================================================================== 359* 360* .. Parameters .. 361 DOUBLE PRECISION ONE 362 PARAMETER ( ONE = 1.0D+0 ) 363* .. 364* .. Local Scalars .. 365 INTEGER I, IINFO, J, JB, NB 366* .. 367* .. External Subroutines .. 368 EXTERNAL DGEMM, DGETF2, DLASWP, DTRSM, XERBLA 369* .. 370* .. External Functions .. 371 INTEGER ILAENV 372 EXTERNAL ILAENV 373* .. 374* .. Intrinsic Functions .. 375 INTRINSIC MAX, MIN 376* .. 377* .. Executable Statements .. 378* 379* Test the input parameters. 380* 381 INFO = 0 382 IF( M.LT.0 ) THEN 383 INFO = -1 384 ELSE IF( N.LT.0 ) THEN 385 INFO = -2 386 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 387 INFO = -4 388 END IF 389 IF( INFO.NE.0 ) THEN 390 CALL XERBLA( 'DGETRF', -INFO ) 391 RETURN 392 END IF 393* 394* Quick return if possible 395* 396 IF( M.EQ.0 .OR. N.EQ.0 ) 397 $ RETURN 398* 399* Determine the block size for this environment. 400* 401 NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 ) 402 IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN 403* 404* Use unblocked code. 405* 406 CALL DGETF2( M, N, A, LDA, IPIV, INFO ) 407 ELSE 408* 409* Use blocked code. 410* 411 DO 20 J = 1, MIN( M, N ), NB 412 JB = MIN( MIN( M, N )-J+1, NB ) 413* 414* Factor diagonal and subdiagonal blocks and test for exact 415* singularity. 416* 417 CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO ) 418* 419* Adjust INFO and the pivot indices. 420* 421 IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 422 $ INFO = IINFO + J - 1 423 DO 10 I = J, MIN( M, J+JB-1 ) 424 IPIV( I ) = J - 1 + IPIV( I ) 425 10 CONTINUE 426* 427* Apply interchanges to columns 1:J-1. 428* 429 CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) 430* 431 IF( J+JB.LE.N ) THEN 432* 433* Apply interchanges to columns J+JB:N. 434* 435 CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, 436 $ IPIV, 1 ) 437* 438* Compute block row of U. 439* 440 CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, 441 $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), 442 $ LDA ) 443 IF( J+JB.LE.M ) THEN 444* 445* Update trailing submatrix. 446* 447 CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1, 448 $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, 449 $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), 450 $ LDA ) 451 END IF 452 END IF 453 20 CONTINUE 454 END IF 455 RETURN 456* 457* End of DGETRF 458* 459 END 460 SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) 461* 462* -- LAPACK routine (version 3.0) -- 463* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 464* Courant Institute, Argonne National Lab, and Rice University 465* March 31, 1993 466* 467* .. Scalar Arguments .. 468 CHARACTER TRANS 469 INTEGER INFO, LDA, LDB, N, NRHS 470* .. 471* .. Array Arguments .. 472 INTEGER IPIV( * ) 473 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 474* .. 475* 476* Purpose 477* ======= 478* 479* DGETRS solves a system of linear equations 480* A * X = B or A' * X = B 481* with a general N-by-N matrix A using the LU factorization computed 482* by DGETRF. 483* 484* Arguments 485* ========= 486* 487* TRANS (input) CHARACTER*1 488* Specifies the form of the system of equations: 489* = 'N': A * X = B (No transpose) 490* = 'T': A'* X = B (Transpose) 491* = 'C': A'* X = B (Conjugate transpose = Transpose) 492* 493* N (input) INTEGER 494* The order of the matrix A. N >= 0. 495* 496* NRHS (input) INTEGER 497* The number of right hand sides, i.e., the number of columns 498* of the matrix B. NRHS >= 0. 499* 500* A (input) DOUBLE PRECISION array, dimension (LDA,N) 501* The factors L and U from the factorization A = P*L*U 502* as computed by DGETRF. 503* 504* LDA (input) INTEGER 505* The leading dimension of the array A. LDA >= max(1,N). 506* 507* IPIV (input) INTEGER array, dimension (N) 508* The pivot indices from DGETRF; for 1<=i<=N, row i of the 509* matrix was interchanged with row IPIV(i). 510* 511* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 512* On entry, the right hand side matrix B. 513* On exit, the solution matrix X. 514* 515* LDB (input) INTEGER 516* The leading dimension of the array B. LDB >= max(1,N). 517* 518* INFO (output) INTEGER 519* = 0: successful exit 520* < 0: if INFO = -i, the i-th argument had an illegal value 521* 522* ===================================================================== 523* 524* .. Parameters .. 525 DOUBLE PRECISION ONE 526 PARAMETER ( ONE = 1.0D+0 ) 527* .. 528* .. Local Scalars .. 529 LOGICAL NOTRAN 530* .. 531* .. External Functions .. 532 LOGICAL LSAME 533 EXTERNAL LSAME 534* .. 535* .. External Subroutines .. 536 EXTERNAL DLASWP, DTRSM, XERBLA 537* .. 538* .. Intrinsic Functions .. 539 INTRINSIC MAX 540* .. 541* .. Executable Statements .. 542* 543* Test the input parameters. 544* 545 INFO = 0 546 NOTRAN = LSAME( TRANS, 'N' ) 547 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 548 $ LSAME( TRANS, 'C' ) ) THEN 549 INFO = -1 550 ELSE IF( N.LT.0 ) THEN 551 INFO = -2 552 ELSE IF( NRHS.LT.0 ) THEN 553 INFO = -3 554 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 555 INFO = -5 556 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 557 INFO = -8 558 END IF 559 IF( INFO.NE.0 ) THEN 560 CALL XERBLA( 'DGETRS', -INFO ) 561 RETURN 562 END IF 563* 564* Quick return if possible 565* 566 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 567 $ RETURN 568* 569 IF( NOTRAN ) THEN 570* 571* Solve A * X = B. 572* 573* Apply row interchanges to the right hand sides. 574* 575 CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, 1 ) 576* 577* Solve L*X = B, overwriting B with X. 578* 579 CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS, 580 $ ONE, A, LDA, B, LDB ) 581* 582* Solve U*X = B, overwriting B with X. 583* 584 CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, 585 $ NRHS, ONE, A, LDA, B, LDB ) 586 ELSE 587* 588* Solve A' * X = B. 589* 590* Solve U'*X = B, overwriting B with X. 591* 592 CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, 593 $ ONE, A, LDA, B, LDB ) 594* 595* Solve L'*X = B, overwriting B with X. 596* 597 CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS, ONE, 598 $ A, LDA, B, LDB ) 599* 600* Apply row interchanges to the solution vectors. 601* 602 CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, -1 ) 603 END IF 604* 605 RETURN 606* 607* End of DGETRS 608* 609 END 610 SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX ) 611* 612* -- LAPACK auxiliary routine (version 3.0) -- 613* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 614* Courant Institute, Argonne National Lab, and Rice University 615* June 30, 1999 616* 617* .. Scalar Arguments .. 618 INTEGER INCX, K1, K2, LDA, N 619* .. 620* .. Array Arguments .. 621 INTEGER IPIV( * ) 622 DOUBLE PRECISION A( LDA, * ) 623* .. 624* 625* Purpose 626* ======= 627* 628* DLASWP performs a series of row interchanges on the matrix A. 629* One row interchange is initiated for each of rows K1 through K2 of A. 630* 631* Arguments 632* ========= 633* 634* N (input) INTEGER 635* The number of columns of the matrix A. 636* 637* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 638* On entry, the matrix of column dimension N to which the row 639* interchanges will be applied. 640* On exit, the permuted matrix. 641* 642* LDA (input) INTEGER 643* The leading dimension of the array A. 644* 645* K1 (input) INTEGER 646* The first element of IPIV for which a row interchange will 647* be done. 648* 649* K2 (input) INTEGER 650* The last element of IPIV for which a row interchange will 651* be done. 652* 653* IPIV (input) INTEGER array, dimension (M*abs(INCX)) 654* The vector of pivot indices. Only the elements in positions 655* K1 through K2 of IPIV are accessed. 656* IPIV(K) = L implies rows K and L are to be interchanged. 657* 658* INCX (input) INTEGER 659* The increment between successive values of IPIV. If IPIV 660* is negative, the pivots are applied in reverse order. 661* 662* Further Details 663* =============== 664* 665* Modified by 666* R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA 667* 668* ===================================================================== 669* 670* .. Local Scalars .. 671 INTEGER I, I1, I2, INC, IP, IX, IX0, J, K, N32 672 DOUBLE PRECISION TEMP 673* .. 674* .. Executable Statements .. 675* 676* Interchange row I with row IPIV(I) for each of rows K1 through K2. 677* 678 IF( INCX.GT.0 ) THEN 679 IX0 = K1 680 I1 = K1 681 I2 = K2 682 INC = 1 683 ELSE IF( INCX.LT.0 ) THEN 684 IX0 = 1 + ( 1-K2 )*INCX 685 I1 = K2 686 I2 = K1 687 INC = -1 688 ELSE 689 RETURN 690 END IF 691* 692 N32 = ( N / 32 )*32 693 IF( N32.NE.0 ) THEN 694 DO 30 J = 1, N32, 32 695 IX = IX0 696 DO 20 I = I1, I2, INC 697 IP = IPIV( IX ) 698 IF( IP.NE.I ) THEN 699 DO 10 K = J, J + 31 700 TEMP = A( I, K ) 701 A( I, K ) = A( IP, K ) 702 A( IP, K ) = TEMP 703 10 CONTINUE 704 END IF 705 IX = IX + INCX 706 20 CONTINUE 707 30 CONTINUE 708 END IF 709 IF( N32.NE.N ) THEN 710 N32 = N32 + 1 711 IX = IX0 712 DO 50 I = I1, I2, INC 713 IP = IPIV( IX ) 714 IF( IP.NE.I ) THEN 715 DO 40 K = N32, N 716 TEMP = A( I, K ) 717 A( I, K ) = A( IP, K ) 718 A( IP, K ) = TEMP 719 40 CONTINUE 720 END IF 721 IX = IX + INCX 722 50 CONTINUE 723 END IF 724* 725 RETURN 726* 727* End of DLASWP 728* 729 END 730 INTEGER FUNCTION IEEECK( ISPEC, ZERO, ONE ) 731* 732* -- LAPACK auxiliary routine (version 3.0) -- 733* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 734* Courant Institute, Argonne National Lab, and Rice University 735* June 30, 1998 736* 737* .. Scalar Arguments .. 738 INTEGER ISPEC 739 REAL ONE, ZERO 740* .. 741* 742* Purpose 743* ======= 744* 745* IEEECK is called from the ILAENV to verify that Infinity and 746* possibly NaN arithmetic is safe (i.e. will not trap). 747* 748* Arguments 749* ========= 750* 751* ISPEC (input) INTEGER 752* Specifies whether to test just for inifinity arithmetic 753* or whether to test for infinity and NaN arithmetic. 754* = 0: Verify infinity arithmetic only. 755* = 1: Verify infinity and NaN arithmetic. 756* 757* ZERO (input) REAL 758* Must contain the value 0.0 759* This is passed to prevent the compiler from optimizing 760* away this code. 761* 762* ONE (input) REAL 763* Must contain the value 1.0 764* This is passed to prevent the compiler from optimizing 765* away this code. 766* 767* RETURN VALUE: INTEGER 768* = 0: Arithmetic failed to produce the correct answers 769* = 1: Arithmetic produced the correct answers 770* 771* .. Local Scalars .. 772 REAL NAN1, NAN2, NAN3, NAN4, NAN5, NAN6, NEGINF, 773 $ NEGZRO, NEWZRO, POSINF 774* .. 775* .. Executable Statements .. 776 IEEECK = 1 777* 778 POSINF = ONE / ZERO 779 IF( POSINF.LE.ONE ) THEN 780 IEEECK = 0 781 RETURN 782 END IF 783* 784 NEGINF = -ONE / ZERO 785 IF( NEGINF.GE.ZERO ) THEN 786 IEEECK = 0 787 RETURN 788 END IF 789* 790 NEGZRO = ONE / ( NEGINF+ONE ) 791 IF( NEGZRO.NE.ZERO ) THEN 792 IEEECK = 0 793 RETURN 794 END IF 795* 796 NEGINF = ONE / NEGZRO 797 IF( NEGINF.GE.ZERO ) THEN 798 IEEECK = 0 799 RETURN 800 END IF 801* 802 NEWZRO = NEGZRO + ZERO 803 IF( NEWZRO.NE.ZERO ) THEN 804 IEEECK = 0 805 RETURN 806 END IF 807* 808 POSINF = ONE / NEWZRO 809 IF( POSINF.LE.ONE ) THEN 810 IEEECK = 0 811 RETURN 812 END IF 813* 814 NEGINF = NEGINF*POSINF 815 IF( NEGINF.GE.ZERO ) THEN 816 IEEECK = 0 817 RETURN 818 END IF 819* 820 POSINF = POSINF*POSINF 821 IF( POSINF.LE.ONE ) THEN 822 IEEECK = 0 823 RETURN 824 END IF 825* 826* 827* 828* 829* Return if we were only asked to check infinity arithmetic 830* 831 IF( ISPEC.EQ.0 ) 832 $ RETURN 833* 834 NAN1 = POSINF + NEGINF 835* 836 NAN2 = POSINF / NEGINF 837* 838 NAN3 = POSINF / POSINF 839* 840 NAN4 = POSINF*ZERO 841* 842 NAN5 = NEGINF*NEGZRO 843* 844 NAN6 = NAN5*0.0 845* 846 IF( NAN1.EQ.NAN1 ) THEN 847 IEEECK = 0 848 RETURN 849 END IF 850* 851 IF( NAN2.EQ.NAN2 ) THEN 852 IEEECK = 0 853 RETURN 854 END IF 855* 856 IF( NAN3.EQ.NAN3 ) THEN 857 IEEECK = 0 858 RETURN 859 END IF 860* 861 IF( NAN4.EQ.NAN4 ) THEN 862 IEEECK = 0 863 RETURN 864 END IF 865* 866 IF( NAN5.EQ.NAN5 ) THEN 867 IEEECK = 0 868 RETURN 869 END IF 870* 871 IF( NAN6.EQ.NAN6 ) THEN 872 IEEECK = 0 873 RETURN 874 END IF 875* 876 RETURN 877 END 878 INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, 879 $ N4 ) 880* 881* -- LAPACK auxiliary routine (version 3.0) -- 882* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 883* Courant Institute, Argonne National Lab, and Rice University 884* June 30, 1999 885* 886* .. Scalar Arguments .. 887 CHARACTER*( * ) NAME, OPTS 888 INTEGER ISPEC, N1, N2, N3, N4 889* .. 890* 891* Purpose 892* ======= 893* 894* ILAENV is called from the LAPACK routines to choose problem-dependent 895* parameters for the local environment. See ISPEC for a description of 896* the parameters. 897* 898* This version provides a set of parameters which should give good, 899* but not optimal, performance on many of the currently available 900* computers. Users are encouraged to modify this subroutine to set 901* the tuning parameters for their particular machine using the option 902* and problem size information in the arguments. 903* 904* This routine will not function correctly if it is converted to all 905* lower case. Converting it to all upper case is allowed. 906* 907* Arguments 908* ========= 909* 910* ISPEC (input) INTEGER 911* Specifies the parameter to be returned as the value of 912* ILAENV. 913* = 1: the optimal blocksize; if this value is 1, an unblocked 914* algorithm will give the best performance. 915* = 2: the minimum block size for which the block routine 916* should be used; if the usable block size is less than 917* this value, an unblocked routine should be used. 918* = 3: the crossover point (in a block routine, for N less 919* than this value, an unblocked routine should be used) 920* = 4: the number of shifts, used in the nonsymmetric 921* eigenvalue routines 922* = 5: the minimum column dimension for blocking to be used; 923* rectangular blocks must have dimension at least k by m, 924* where k is given by ILAENV(2,...) and m by ILAENV(5,...) 925* = 6: the crossover point for the SVD (when reducing an m by n 926* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds 927* this value, a QR factorization is used first to reduce 928* the matrix to a triangular form.) 929* = 7: the number of processors 930* = 8: the crossover point for the multishift QR and QZ methods 931* for nonsymmetric eigenvalue problems. 932* = 9: maximum size of the subproblems at the bottom of the 933* computation tree in the divide-and-conquer algorithm 934* (used by xGELSD and xGESDD) 935* =10: ieee NaN arithmetic can be trusted not to trap 936* =11: infinity arithmetic can be trusted not to trap 937* 938* NAME (input) CHARACTER*(*) 939* The name of the calling subroutine, in either upper case or 940* lower case. 941* 942* OPTS (input) CHARACTER*(*) 943* The character options to the subroutine NAME, concatenated 944* into a single character string. For example, UPLO = 'U', 945* TRANS = 'T', and DIAG = 'N' for a triangular routine would 946* be specified as OPTS = 'UTN'. 947* 948* N1 (input) INTEGER 949* N2 (input) INTEGER 950* N3 (input) INTEGER 951* N4 (input) INTEGER 952* Problem dimensions for the subroutine NAME; these may not all 953* be required. 954* 955* (ILAENV) (output) INTEGER 956* >= 0: the value of the parameter specified by ISPEC 957* < 0: if ILAENV = -k, the k-th argument had an illegal value. 958* 959* Further Details 960* =============== 961* 962* The following conventions have been used when calling ILAENV from the 963* LAPACK routines: 964* 1) OPTS is a concatenation of all of the character options to 965* subroutine NAME, in the same order that they appear in the 966* argument list for NAME, even if they are not used in determining 967* the value of the parameter specified by ISPEC. 968* 2) The problem dimensions N1, N2, N3, N4 are specified in the order 969* that they appear in the argument list for NAME. N1 is used 970* first, N2 second, and so on, and unused problem dimensions are 971* passed a value of -1. 972* 3) The parameter value returned by ILAENV is checked for validity in 973* the calling subroutine. For example, ILAENV is used to retrieve 974* the optimal blocksize for STRTRI as follows: 975* 976* NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) 977* IF( NB.LE.1 ) NB = MAX( 1, N ) 978* 979* ===================================================================== 980* 981* .. Local Scalars .. 982 LOGICAL CNAME, SNAME 983 CHARACTER*1 C1 984 CHARACTER*2 C2, C4 985 CHARACTER*3 C3 986 CHARACTER*6 SUBNAM 987 INTEGER I, IC, IZ, NB, NBMIN, NX 988* .. 989* .. Intrinsic Functions .. 990 INTRINSIC CHAR, ICHAR, INT, MIN, REAL 991* .. 992* .. External Functions .. 993 INTEGER IEEECK 994 EXTERNAL IEEECK 995* .. 996* .. Executable Statements .. 997* 998C GO TO ( 100, 100, 100, 400, 500, 600, 700, 800, 900, 1000, 999C $ 1100 ) ISPEC 1000C CHANGED COMPUTED GOTO: OBSOLETE 1001C 1002 IF((ISPEC.EQ.1).OR.(ISPEC.EQ.2).OR.(ISPEC.EQ.3)) THEN 1003 GO TO 100 1004 ELSEIF(ISPEC.EQ.4) THEN 1005 GO TO 400 1006 ELSEIF(ISPEC.EQ.5) THEN 1007 GO TO 500 1008 ELSEIF(ISPEC.EQ.6) THEN 1009 GO TO 600 1010 ELSEIF(ISPEC.EQ.7) THEN 1011 GO TO 700 1012 ELSEIF(ISPEC.EQ.8) THEN 1013 GO TO 800 1014 ELSEIF(ISPEC.EQ.9) THEN 1015 GO TO 900 1016 ELSEIF(ISPEC.EQ.10) THEN 1017 GO TO 1000 1018 ELSEIF(ISPEC.EQ.11) THEN 1019 GO TO 1100 1020 ENDIF 1021* 1022* Invalid value for ISPEC 1023* 1024 ILAENV = -1 1025 RETURN 1026* 1027 100 CONTINUE 1028* 1029* Convert NAME to upper case if the first character is lower case. 1030* 1031 ILAENV = 1 1032 SUBNAM = NAME 1033 IC = ICHAR( SUBNAM( 1:1 ) ) 1034 IZ = ICHAR( 'Z' ) 1035 IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN 1036* 1037* ASCII character set 1038* 1039 IF( IC.GE.97 .AND. IC.LE.122 ) THEN 1040 SUBNAM( 1:1 ) = CHAR( IC-32 ) 1041 DO 10 I = 2, 6 1042 IC = ICHAR( SUBNAM( I:I ) ) 1043 IF( IC.GE.97 .AND. IC.LE.122 ) 1044 $ SUBNAM( I:I ) = CHAR( IC-32 ) 1045 10 CONTINUE 1046 END IF 1047* 1048 ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN 1049* 1050* EBCDIC character set 1051* 1052 IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. 1053 $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. 1054 $ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN 1055 SUBNAM( 1:1 ) = CHAR( IC+64 ) 1056 DO 20 I = 2, 6 1057 IC = ICHAR( SUBNAM( I:I ) ) 1058 IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. 1059 $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. 1060 $ ( IC.GE.162 .AND. IC.LE.169 ) ) 1061 $ SUBNAM( I:I ) = CHAR( IC+64 ) 1062 20 CONTINUE 1063 END IF 1064* 1065 ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN 1066* 1067* Prime machines: ASCII+128 1068* 1069 IF( IC.GE.225 .AND. IC.LE.250 ) THEN 1070 SUBNAM( 1:1 ) = CHAR( IC-32 ) 1071 DO 30 I = 2, 6 1072 IC = ICHAR( SUBNAM( I:I ) ) 1073 IF( IC.GE.225 .AND. IC.LE.250 ) 1074 $ SUBNAM( I:I ) = CHAR( IC-32 ) 1075 30 CONTINUE 1076 END IF 1077 END IF 1078* 1079 C1 = SUBNAM( 1:1 ) 1080 SNAME = C1.EQ.'S' .OR. C1.EQ.'D' 1081 CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' 1082 IF( .NOT.( CNAME .OR. SNAME ) ) 1083 $ RETURN 1084 C2 = SUBNAM( 2:3 ) 1085 C3 = SUBNAM( 4:6 ) 1086 C4 = C3( 2:3 ) 1087* 1088 GO TO ( 110, 200, 300 ) ISPEC 1089* 1090 110 CONTINUE 1091* 1092* ISPEC = 1: block size 1093* 1094* In these examples, separate code is provided for setting NB for 1095* real and complex. We assume that NB will take the same value in 1096* single or double precision. 1097* 1098 NB = 1 1099* 1100 IF( C2.EQ.'GE' ) THEN 1101 IF( C3.EQ.'TRF' ) THEN 1102 IF( SNAME ) THEN 1103 NB = 64 1104 ELSE 1105 NB = 64 1106 END IF 1107 ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. 1108 $ C3.EQ.'QLF' ) THEN 1109 IF( SNAME ) THEN 1110 NB = 32 1111 ELSE 1112 NB = 32 1113 END IF 1114 ELSE IF( C3.EQ.'HRD' ) THEN 1115 IF( SNAME ) THEN 1116 NB = 32 1117 ELSE 1118 NB = 32 1119 END IF 1120 ELSE IF( C3.EQ.'BRD' ) THEN 1121 IF( SNAME ) THEN 1122 NB = 32 1123 ELSE 1124 NB = 32 1125 END IF 1126 ELSE IF( C3.EQ.'TRI' ) THEN 1127 IF( SNAME ) THEN 1128 NB = 64 1129 ELSE 1130 NB = 64 1131 END IF 1132 END IF 1133 ELSE IF( C2.EQ.'PO' ) THEN 1134 IF( C3.EQ.'TRF' ) THEN 1135 IF( SNAME ) THEN 1136 NB = 64 1137 ELSE 1138 NB = 64 1139 END IF 1140 END IF 1141 ELSE IF( C2.EQ.'SY' ) THEN 1142 IF( C3.EQ.'TRF' ) THEN 1143 IF( SNAME ) THEN 1144 NB = 64 1145 ELSE 1146 NB = 64 1147 END IF 1148 ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN 1149 NB = 32 1150 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN 1151 NB = 64 1152 END IF 1153 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 1154 IF( C3.EQ.'TRF' ) THEN 1155 NB = 64 1156 ELSE IF( C3.EQ.'TRD' ) THEN 1157 NB = 32 1158 ELSE IF( C3.EQ.'GST' ) THEN 1159 NB = 64 1160 END IF 1161 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 1162 IF( C3( 1:1 ).EQ.'G' ) THEN 1163 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1164 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1165 $ C4.EQ.'BR' ) THEN 1166 NB = 32 1167 END IF 1168 ELSE IF( C3( 1:1 ).EQ.'M' ) THEN 1169 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1170 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1171 $ C4.EQ.'BR' ) THEN 1172 NB = 32 1173 END IF 1174 END IF 1175 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 1176 IF( C3( 1:1 ).EQ.'G' ) THEN 1177 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1178 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1179 $ C4.EQ.'BR' ) THEN 1180 NB = 32 1181 END IF 1182 ELSE IF( C3( 1:1 ).EQ.'M' ) THEN 1183 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1184 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1185 $ C4.EQ.'BR' ) THEN 1186 NB = 32 1187 END IF 1188 END IF 1189 ELSE IF( C2.EQ.'GB' ) THEN 1190 IF( C3.EQ.'TRF' ) THEN 1191 IF( SNAME ) THEN 1192 IF( N4.LE.64 ) THEN 1193 NB = 1 1194 ELSE 1195 NB = 32 1196 END IF 1197 ELSE 1198 IF( N4.LE.64 ) THEN 1199 NB = 1 1200 ELSE 1201 NB = 32 1202 END IF 1203 END IF 1204 END IF 1205 ELSE IF( C2.EQ.'PB' ) THEN 1206 IF( C3.EQ.'TRF' ) THEN 1207 IF( SNAME ) THEN 1208 IF( N2.LE.64 ) THEN 1209 NB = 1 1210 ELSE 1211 NB = 32 1212 END IF 1213 ELSE 1214 IF( N2.LE.64 ) THEN 1215 NB = 1 1216 ELSE 1217 NB = 32 1218 END IF 1219 END IF 1220 END IF 1221 ELSE IF( C2.EQ.'TR' ) THEN 1222 IF( C3.EQ.'TRI' ) THEN 1223 IF( SNAME ) THEN 1224 NB = 64 1225 ELSE 1226 NB = 64 1227 END IF 1228 END IF 1229 ELSE IF( C2.EQ.'LA' ) THEN 1230 IF( C3.EQ.'UUM' ) THEN 1231 IF( SNAME ) THEN 1232 NB = 64 1233 ELSE 1234 NB = 64 1235 END IF 1236 END IF 1237 ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN 1238 IF( C3.EQ.'EBZ' ) THEN 1239 NB = 1 1240 END IF 1241 END IF 1242 ILAENV = NB 1243 RETURN 1244* 1245 200 CONTINUE 1246* 1247* ISPEC = 2: minimum block size 1248* 1249 NBMIN = 2 1250 IF( C2.EQ.'GE' ) THEN 1251 IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. 1252 $ C3.EQ.'QLF' ) THEN 1253 IF( SNAME ) THEN 1254 NBMIN = 2 1255 ELSE 1256 NBMIN = 2 1257 END IF 1258 ELSE IF( C3.EQ.'HRD' ) THEN 1259 IF( SNAME ) THEN 1260 NBMIN = 2 1261 ELSE 1262 NBMIN = 2 1263 END IF 1264 ELSE IF( C3.EQ.'BRD' ) THEN 1265 IF( SNAME ) THEN 1266 NBMIN = 2 1267 ELSE 1268 NBMIN = 2 1269 END IF 1270 ELSE IF( C3.EQ.'TRI' ) THEN 1271 IF( SNAME ) THEN 1272 NBMIN = 2 1273 ELSE 1274 NBMIN = 2 1275 END IF 1276 END IF 1277 ELSE IF( C2.EQ.'SY' ) THEN 1278 IF( C3.EQ.'TRF' ) THEN 1279 IF( SNAME ) THEN 1280 NBMIN = 8 1281 ELSE 1282 NBMIN = 8 1283 END IF 1284 ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN 1285 NBMIN = 2 1286 END IF 1287 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 1288 IF( C3.EQ.'TRD' ) THEN 1289 NBMIN = 2 1290 END IF 1291 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 1292 IF( C3( 1:1 ).EQ.'G' ) THEN 1293 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1294 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1295 $ C4.EQ.'BR' ) THEN 1296 NBMIN = 2 1297 END IF 1298 ELSE IF( C3( 1:1 ).EQ.'M' ) THEN 1299 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1300 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1301 $ C4.EQ.'BR' ) THEN 1302 NBMIN = 2 1303 END IF 1304 END IF 1305 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 1306 IF( C3( 1:1 ).EQ.'G' ) THEN 1307 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1308 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1309 $ C4.EQ.'BR' ) THEN 1310 NBMIN = 2 1311 END IF 1312 ELSE IF( C3( 1:1 ).EQ.'M' ) THEN 1313 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1314 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1315 $ C4.EQ.'BR' ) THEN 1316 NBMIN = 2 1317 END IF 1318 END IF 1319 END IF 1320 ILAENV = NBMIN 1321 RETURN 1322* 1323 300 CONTINUE 1324* 1325* ISPEC = 3: crossover point 1326* 1327 NX = 0 1328 IF( C2.EQ.'GE' ) THEN 1329 IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. 1330 $ C3.EQ.'QLF' ) THEN 1331 IF( SNAME ) THEN 1332 NX = 128 1333 ELSE 1334 NX = 128 1335 END IF 1336 ELSE IF( C3.EQ.'HRD' ) THEN 1337 IF( SNAME ) THEN 1338 NX = 128 1339 ELSE 1340 NX = 128 1341 END IF 1342 ELSE IF( C3.EQ.'BRD' ) THEN 1343 IF( SNAME ) THEN 1344 NX = 128 1345 ELSE 1346 NX = 128 1347 END IF 1348 END IF 1349 ELSE IF( C2.EQ.'SY' ) THEN 1350 IF( SNAME .AND. C3.EQ.'TRD' ) THEN 1351 NX = 32 1352 END IF 1353 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 1354 IF( C3.EQ.'TRD' ) THEN 1355 NX = 32 1356 END IF 1357 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 1358 IF( C3( 1:1 ).EQ.'G' ) THEN 1359 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1360 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1361 $ C4.EQ.'BR' ) THEN 1362 NX = 128 1363 END IF 1364 END IF 1365 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 1366 IF( C3( 1:1 ).EQ.'G' ) THEN 1367 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 1368 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 1369 $ C4.EQ.'BR' ) THEN 1370 NX = 128 1371 END IF 1372 END IF 1373 END IF 1374 ILAENV = NX 1375 RETURN 1376* 1377 400 CONTINUE 1378* 1379* ISPEC = 4: number of shifts (used by xHSEQR) 1380* 1381 ILAENV = 6 1382 RETURN 1383* 1384 500 CONTINUE 1385* 1386* ISPEC = 5: minimum column dimension (not used) 1387* 1388 ILAENV = 2 1389 RETURN 1390* 1391 600 CONTINUE 1392* 1393* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) 1394* 1395 ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) 1396 RETURN 1397* 1398 700 CONTINUE 1399* 1400* ISPEC = 7: number of processors (not used) 1401* 1402 ILAENV = 1 1403 RETURN 1404* 1405 800 CONTINUE 1406* 1407* ISPEC = 8: crossover point for multishift (used by xHSEQR) 1408* 1409 ILAENV = 50 1410 RETURN 1411* 1412 900 CONTINUE 1413* 1414* ISPEC = 9: maximum size of the subproblems at the bottom of the 1415* computation tree in the divide-and-conquer algorithm 1416* (used by xGELSD and xGESDD) 1417* 1418 ILAENV = 25 1419 RETURN 1420* 1421 1000 CONTINUE 1422* 1423* ISPEC = 10: ieee NaN arithmetic can be trusted not to trap 1424* 1425c ILAENV = 0 1426 ILAENV = 1 1427 IF( ILAENV.EQ.1 ) THEN 1428 ILAENV = IEEECK( 0, 0.0, 1.0 ) 1429 END IF 1430 RETURN 1431* 1432 1100 CONTINUE 1433* 1434* ISPEC = 11: infinity arithmetic can be trusted not to trap 1435* 1436c ILAENV = 0 1437 ILAENV = 1 1438 IF( ILAENV.EQ.1 ) THEN 1439 ILAENV = IEEECK( 1, 0.0, 1.0 ) 1440 END IF 1441 RETURN 1442* 1443* End of ILAENV 1444* 1445 END 1446cc 1447 integer function idamax(n,dx,incx) 1448c 1449c finds the index of element having max. absolute value. 1450c jack dongarra, linpack, 3/11/78. 1451c modified 3/93 to return if incx .le. 0. 1452c modified 12/3/93, array(1) declarations changed to array(*) 1453c 1454 double precision dx(*),dmax 1455 integer i,incx,ix,n 1456c 1457 idamax = 0 1458 if( n.lt.1 .or. incx.le.0 ) return 1459 idamax = 1 1460 if(n.eq.1)return 1461 if(incx.eq.1)go to 20 1462c 1463c code for increment not equal to 1 1464c 1465 ix = 1 1466 dmax = dabs(dx(1)) 1467 ix = ix + incx 1468 do 10 i = 2,n 1469 if(dabs(dx(ix)).le.dmax) go to 5 1470 idamax = i 1471 dmax = dabs(dx(ix)) 1472 5 ix = ix + incx 1473 10 continue 1474 return 1475c 1476c code for increment equal to 1 1477c 1478 20 dmax = dabs(dx(1)) 1479 do 30 i = 2,n 1480 if(dabs(dx(i)).le.dmax) go to 30 1481 idamax = i 1482 dmax = dabs(dx(i)) 1483 30 continue 1484 return 1485 end 1486cc 1487 SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) 1488* .. Scalar Arguments .. 1489 DOUBLE PRECISION ALPHA 1490 INTEGER INCX, INCY, LDA, M, N 1491* .. Array Arguments .. 1492 DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) 1493* .. 1494* 1495* Purpose 1496* ======= 1497* 1498* DGER performs the rank 1 operation 1499* 1500* A := alpha*x*y' + A, 1501* 1502* where alpha is a scalar, x is an m element vector, y is an n element 1503* vector and A is an m by n matrix. 1504* 1505* Parameters 1506* ========== 1507* 1508* M - INTEGER. 1509* On entry, M specifies the number of rows of the matrix A. 1510* M must be at least zero. 1511* Unchanged on exit. 1512* 1513* N - INTEGER. 1514* On entry, N specifies the number of columns of the matrix A. 1515* N must be at least zero. 1516* Unchanged on exit. 1517* 1518* ALPHA - DOUBLE PRECISION. 1519* On entry, ALPHA specifies the scalar alpha. 1520* Unchanged on exit. 1521* 1522* X - DOUBLE PRECISION array of dimension at least 1523* ( 1 + ( m - 1 )*abs( INCX ) ). 1524* Before entry, the incremented array X must contain the m 1525* element vector x. 1526* Unchanged on exit. 1527* 1528* INCX - INTEGER. 1529* On entry, INCX specifies the increment for the elements of 1530* X. INCX must not be zero. 1531* Unchanged on exit. 1532* 1533* Y - DOUBLE PRECISION array of dimension at least 1534* ( 1 + ( n - 1 )*abs( INCY ) ). 1535* Before entry, the incremented array Y must contain the n 1536* element vector y. 1537* Unchanged on exit. 1538* 1539* INCY - INTEGER. 1540* On entry, INCY specifies the increment for the elements of 1541* Y. INCY must not be zero. 1542* Unchanged on exit. 1543* 1544* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 1545* Before entry, the leading m by n part of the array A must 1546* contain the matrix of coefficients. On exit, A is 1547* overwritten by the updated matrix. 1548* 1549* LDA - INTEGER. 1550* On entry, LDA specifies the first dimension of A as declared 1551* in the calling (sub) program. LDA must be at least 1552* max( 1, m ). 1553* Unchanged on exit. 1554* 1555* 1556* Level 2 Blas routine. 1557* 1558* -- Written on 22-October-1986. 1559* Jack Dongarra, Argonne National Lab. 1560* Jeremy Du Croz, Nag Central Office. 1561* Sven Hammarling, Nag Central Office. 1562* Richard Hanson, Sandia National Labs. 1563* 1564* 1565* .. Parameters .. 1566 DOUBLE PRECISION ZERO 1567 PARAMETER ( ZERO = 0.0D+0 ) 1568* .. Local Scalars .. 1569 DOUBLE PRECISION TEMP 1570 INTEGER I, INFO, IX, J, JY, KX 1571* .. External Subroutines .. 1572 EXTERNAL XERBLA 1573* .. Intrinsic Functions .. 1574 INTRINSIC MAX 1575* .. 1576* .. Executable Statements .. 1577* 1578* Test the input parameters. 1579* 1580 INFO = 0 1581 IF ( M.LT.0 )THEN 1582 INFO = 1 1583 ELSE IF( N.LT.0 )THEN 1584 INFO = 2 1585 ELSE IF( INCX.EQ.0 )THEN 1586 INFO = 5 1587 ELSE IF( INCY.EQ.0 )THEN 1588 INFO = 7 1589 ELSE IF( LDA.LT.MAX( 1, M ) )THEN 1590 INFO = 9 1591 END IF 1592 IF( INFO.NE.0 )THEN 1593 CALL XERBLA( 'DGER ', INFO ) 1594 RETURN 1595 END IF 1596* 1597* Quick return if possible. 1598* 1599 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) 1600 $ RETURN 1601* 1602* Start the operations. In this version the elements of A are 1603* accessed sequentially with one pass through A. 1604* 1605 IF( INCY.GT.0 )THEN 1606 JY = 1 1607 ELSE 1608 JY = 1 - ( N - 1 )*INCY 1609 END IF 1610 IF( INCX.EQ.1 )THEN 1611 DO 20, J = 1, N 1612 IF( Y( JY ).NE.ZERO )THEN 1613 TEMP = ALPHA*Y( JY ) 1614 DO 10, I = 1, M 1615 A( I, J ) = A( I, J ) + X( I )*TEMP 1616 10 CONTINUE 1617 END IF 1618 JY = JY + INCY 1619 20 CONTINUE 1620 ELSE 1621 IF( INCX.GT.0 )THEN 1622 KX = 1 1623 ELSE 1624 KX = 1 - ( M - 1 )*INCX 1625 END IF 1626 DO 40, J = 1, N 1627 IF( Y( JY ).NE.ZERO )THEN 1628 TEMP = ALPHA*Y( JY ) 1629 IX = KX 1630 DO 30, I = 1, M 1631 A( I, J ) = A( I, J ) + X( IX )*TEMP 1632 IX = IX + INCX 1633 30 CONTINUE 1634 END IF 1635 JY = JY + INCY 1636 40 CONTINUE 1637 END IF 1638* 1639 RETURN 1640* 1641* End of DGER . 1642* 1643 END 1644c SUBROUTINE XERBLA( SRNAME, INFO ) 1645c* 1646c* -- LAPACK auxiliary routine (preliminary version) -- 1647c* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 1648c* Courant Institute, Argonne National Lab, and Rice University 1649c* February 29, 1992 1650c* 1651c* .. Scalar Arguments .. 1652c CHARACTER*6 SRNAME 1653c INTEGER INFO 1654c* .. 1655c* 1656c* Purpose 1657c* ======= 1658c* 1659c* XERBLA is an error handler for the LAPACK routines. 1660c* It is called by an LAPACK routine if an input parameter has an 1661c* invalid value. A message is printed and execution stops. 1662c* 1663c* Installers may consider modifying the STOP statement in order to 1664c* call system-specific exception-handling facilities. 1665cc* 1666c* Arguments 1667c* ========= 1668c* 1669c* SRNAME (input) CHARACTER*6 1670c* The name of the routine which called XERBLA. 1671c* 1672c* INFO (input) INTEGER 1673c* The position of the invalid parameter in the parameter list 1674c* of the calling routine. 1675c* 1676c* 1677c WRITE( *, FMT = 9999 )SRNAME, INFO 1678c* 1679c STOP 1680c* 1681c 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', 1682c $ 'an illegal value' ) 1683c* 1684c* End of XERBLA 1685c* 1686c END 1687cc 1688 subroutine dswap (n,dx,incx,dy,incy) 1689c 1690c interchanges two vectors. 1691c uses unrolled loops for increments equal one. 1692c jack dongarra, linpack, 3/11/78. 1693c modified 12/3/93, array(1) declarations changed to array(*) 1694c 1695 double precision dx(*),dy(*),dtemp 1696 integer i,incx,incy,ix,iy,m,mp1,n 1697c 1698 if(n.le.0)return 1699 if(incx.eq.1.and.incy.eq.1)go to 20 1700c 1701c code for unequal increments or equal increments not equal 1702c to 1 1703c 1704 ix = 1 1705 iy = 1 1706 if(incx.lt.0)ix = (-n+1)*incx + 1 1707 if(incy.lt.0)iy = (-n+1)*incy + 1 1708 do 10 i = 1,n 1709 dtemp = dx(ix) 1710 dx(ix) = dy(iy) 1711 dy(iy) = dtemp 1712 ix = ix + incx 1713 iy = iy + incy 1714 10 continue 1715 return 1716c 1717c code for both increments equal to 1 1718c 1719c 1720c clean-up loop 1721c 1722 20 m = mod(n,3) 1723 if( m .eq. 0 ) go to 40 1724 do 30 i = 1,m 1725 dtemp = dx(i) 1726 dx(i) = dy(i) 1727 dy(i) = dtemp 1728 30 continue 1729 if( n .lt. 3 ) return 1730 40 mp1 = m + 1 1731 do 50 i = mp1,n,3 1732 dtemp = dx(i) 1733 dx(i) = dy(i) 1734 dy(i) = dtemp 1735 dtemp = dx(i + 1) 1736 dx(i + 1) = dy(i + 1) 1737 dy(i + 1) = dtemp 1738 dtemp = dx(i + 2) 1739 dx(i + 2) = dy(i + 2) 1740 dy(i + 2) = dtemp 1741 50 continue 1742 return 1743 end 1744cc 1745 subroutine dscal(n,da,dx,incx) 1746c 1747c scales a vector by a constant. 1748c uses unrolled loops for increment equal to one. 1749c jack dongarra, linpack, 3/11/78. 1750c modified 3/93 to return if incx .le. 0. 1751c modified 12/3/93, array(1) declarations changed to array(*) 1752c 1753 double precision da,dx(*) 1754 integer i,incx,m,mp1,n,nincx 1755c 1756 if( n.le.0 .or. incx.le.0 )return 1757 if(incx.eq.1)go to 20 1758c 1759c code for increment not equal to 1 1760c 1761 nincx = n*incx 1762 do 10 i = 1,nincx,incx 1763 dx(i) = da*dx(i) 1764 10 continue 1765 return 1766c 1767c code for increment equal to 1 1768c 1769c 1770c clean-up loop 1771c 1772 20 m = mod(n,5) 1773 if( m .eq. 0 ) go to 40 1774 do 30 i = 1,m 1775 dx(i) = da*dx(i) 1776 30 continue 1777 if( n .lt. 5 ) return 1778 40 mp1 = m + 1 1779 do 50 i = mp1,n,5 1780 dx(i) = da*dx(i) 1781 dx(i + 1) = da*dx(i + 1) 1782 dx(i + 2) = da*dx(i + 2) 1783 dx(i + 3) = da*dx(i + 3) 1784 dx(i + 4) = da*dx(i + 4) 1785 50 continue 1786 return 1787 end 1788cc 1789 SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, 1790 $ B, LDB ) 1791* .. Scalar Arguments .. 1792 CHARACTER*1 SIDE, UPLO, TRANSA, DIAG 1793 INTEGER M, N, LDA, LDB 1794 DOUBLE PRECISION ALPHA 1795* .. Array Arguments .. 1796 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 1797* .. 1798* 1799* Purpose 1800* ======= 1801* 1802* DTRSM solves one of the matrix equations 1803* 1804* op( A )*X = alpha*B, or X*op( A ) = alpha*B, 1805* 1806* where alpha is a scalar, X and B are m by n matrices, A is a unit, or 1807* non-unit, upper or lower triangular matrix and op( A ) is one of 1808* 1809* op( A ) = A or op( A ) = A'. 1810* 1811* The matrix X is overwritten on B. 1812* 1813* Parameters 1814* ========== 1815* 1816* SIDE - CHARACTER*1. 1817* On entry, SIDE specifies whether op( A ) appears on the left 1818* or right of X as follows: 1819* 1820* SIDE = 'L' or 'l' op( A )*X = alpha*B. 1821* 1822* SIDE = 'R' or 'r' X*op( A ) = alpha*B. 1823* 1824* Unchanged on exit. 1825* 1826* UPLO - CHARACTER*1. 1827* On entry, UPLO specifies whether the matrix A is an upper or 1828* lower triangular matrix as follows: 1829* 1830* UPLO = 'U' or 'u' A is an upper triangular matrix. 1831* 1832* UPLO = 'L' or 'l' A is a lower triangular matrix. 1833* 1834* Unchanged on exit. 1835* 1836* TRANSA - CHARACTER*1. 1837* On entry, TRANSA specifies the form of op( A ) to be used in 1838* the matrix multiplication as follows: 1839* 1840* TRANSA = 'N' or 'n' op( A ) = A. 1841* 1842* TRANSA = 'T' or 't' op( A ) = A'. 1843* 1844* TRANSA = 'C' or 'c' op( A ) = A'. 1845* 1846* Unchanged on exit. 1847* 1848* DIAG - CHARACTER*1. 1849* On entry, DIAG specifies whether or not A is unit triangular 1850* as follows: 1851* 1852* DIAG = 'U' or 'u' A is assumed to be unit triangular. 1853* 1854* DIAG = 'N' or 'n' A is not assumed to be unit 1855* triangular. 1856* 1857* Unchanged on exit. 1858* 1859* M - INTEGER. 1860* On entry, M specifies the number of rows of B. M must be at 1861* least zero. 1862* Unchanged on exit. 1863* 1864* N - INTEGER. 1865* On entry, N specifies the number of columns of B. N must be 1866* at least zero. 1867* Unchanged on exit. 1868* 1869* ALPHA - DOUBLE PRECISION. 1870* On entry, ALPHA specifies the scalar alpha. When alpha is 1871* zero then A is not referenced and B need not be set before 1872* entry. 1873* Unchanged on exit. 1874* 1875* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m 1876* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. 1877* Before entry with UPLO = 'U' or 'u', the leading k by k 1878* upper triangular part of the array A must contain the upper 1879* triangular matrix and the strictly lower triangular part of 1880* A is not referenced. 1881* Before entry with UPLO = 'L' or 'l', the leading k by k 1882* lower triangular part of the array A must contain the lower 1883* triangular matrix and the strictly upper triangular part of 1884* A is not referenced. 1885* Note that when DIAG = 'U' or 'u', the diagonal elements of 1886* A are not referenced either, but are assumed to be unity. 1887* Unchanged on exit. 1888* 1889* LDA - INTEGER. 1890* On entry, LDA specifies the first dimension of A as declared 1891* in the calling (sub) program. When SIDE = 'L' or 'l' then 1892* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' 1893* then LDA must be at least max( 1, n ). 1894* Unchanged on exit. 1895* 1896* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). 1897* Before entry, the leading m by n part of the array B must 1898* contain the right-hand side matrix B, and on exit is 1899* overwritten by the solution matrix X. 1900* 1901* LDB - INTEGER. 1902* On entry, LDB specifies the first dimension of B as declared 1903* in the calling (sub) program. LDB must be at least 1904* max( 1, m ). 1905* Unchanged on exit. 1906* 1907* 1908* Level 3 Blas routine. 1909* 1910* 1911* -- Written on 8-February-1989. 1912* Jack Dongarra, Argonne National Laboratory. 1913* Iain Duff, AERE Harwell. 1914* Jeremy Du Croz, Numerical Algorithms Group Ltd. 1915* Sven Hammarling, Numerical Algorithms Group Ltd. 1916* 1917* 1918* .. External Functions .. 1919 LOGICAL LSAME 1920 EXTERNAL LSAME 1921* .. External Subroutines .. 1922 EXTERNAL XERBLA 1923* .. Intrinsic Functions .. 1924 INTRINSIC MAX 1925* .. Local Scalars .. 1926 LOGICAL LSIDE, NOUNIT, UPPER 1927 INTEGER I, INFO, J, K, NROWA 1928 DOUBLE PRECISION TEMP 1929* .. Parameters .. 1930 DOUBLE PRECISION ONE , ZERO 1931 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 1932* .. 1933* .. Executable Statements .. 1934* 1935* Test the input parameters. 1936* 1937 LSIDE = LSAME( SIDE , 'L' ) 1938 IF( LSIDE )THEN 1939 NROWA = M 1940 ELSE 1941 NROWA = N 1942 END IF 1943 NOUNIT = LSAME( DIAG , 'N' ) 1944 UPPER = LSAME( UPLO , 'U' ) 1945* 1946 INFO = 0 1947 IF( ( .NOT.LSIDE ).AND. 1948 $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN 1949 INFO = 1 1950 ELSE IF( ( .NOT.UPPER ).AND. 1951 $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN 1952 INFO = 2 1953 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. 1954 $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. 1955 $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN 1956 INFO = 3 1957 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. 1958 $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN 1959 INFO = 4 1960 ELSE IF( M .LT.0 )THEN 1961 INFO = 5 1962 ELSE IF( N .LT.0 )THEN 1963 INFO = 6 1964 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN 1965 INFO = 9 1966 ELSE IF( LDB.LT.MAX( 1, M ) )THEN 1967 INFO = 11 1968 END IF 1969 IF( INFO.NE.0 )THEN 1970 CALL XERBLA( 'DTRSM ', INFO ) 1971 RETURN 1972 END IF 1973* 1974* Quick return if possible. 1975* 1976 IF( N.EQ.0 ) 1977 $ RETURN 1978* 1979* And when alpha.eq.zero. 1980* 1981 IF( ALPHA.EQ.ZERO )THEN 1982 DO 20, J = 1, N 1983 DO 10, I = 1, M 1984 B( I, J ) = ZERO 1985 10 CONTINUE 1986 20 CONTINUE 1987 RETURN 1988 END IF 1989* 1990* Start the operations. 1991* 1992 IF( LSIDE )THEN 1993 IF( LSAME( TRANSA, 'N' ) )THEN 1994* 1995* Form B := alpha*inv( A )*B. 1996* 1997 IF( UPPER )THEN 1998 DO 60, J = 1, N 1999 IF( ALPHA.NE.ONE )THEN 2000 DO 30, I = 1, M 2001 B( I, J ) = ALPHA*B( I, J ) 2002 30 CONTINUE 2003 END IF 2004 DO 50, K = M, 1, -1 2005 IF( B( K, J ).NE.ZERO )THEN 2006 IF( NOUNIT ) 2007 $ B( K, J ) = B( K, J )/A( K, K ) 2008 DO 40, I = 1, K - 1 2009 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 2010 40 CONTINUE 2011 END IF 2012 50 CONTINUE 2013 60 CONTINUE 2014 ELSE 2015 DO 100, J = 1, N 2016 IF( ALPHA.NE.ONE )THEN 2017 DO 70, I = 1, M 2018 B( I, J ) = ALPHA*B( I, J ) 2019 70 CONTINUE 2020 END IF 2021 DO 90 K = 1, M 2022 IF( B( K, J ).NE.ZERO )THEN 2023 IF( NOUNIT ) 2024 $ B( K, J ) = B( K, J )/A( K, K ) 2025 DO 80, I = K + 1, M 2026 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 2027 80 CONTINUE 2028 END IF 2029 90 CONTINUE 2030 100 CONTINUE 2031 END IF 2032 ELSE 2033* 2034* Form B := alpha*inv( A' )*B. 2035* 2036 IF( UPPER )THEN 2037 DO 130, J = 1, N 2038 DO 120, I = 1, M 2039 TEMP = ALPHA*B( I, J ) 2040 DO 110, K = 1, I - 1 2041 TEMP = TEMP - A( K, I )*B( K, J ) 2042 110 CONTINUE 2043 IF( NOUNIT ) 2044 $ TEMP = TEMP/A( I, I ) 2045 B( I, J ) = TEMP 2046 120 CONTINUE 2047 130 CONTINUE 2048 ELSE 2049 DO 160, J = 1, N 2050 DO 150, I = M, 1, -1 2051 TEMP = ALPHA*B( I, J ) 2052 DO 140, K = I + 1, M 2053 TEMP = TEMP - A( K, I )*B( K, J ) 2054 140 CONTINUE 2055 IF( NOUNIT ) 2056 $ TEMP = TEMP/A( I, I ) 2057 B( I, J ) = TEMP 2058 150 CONTINUE 2059 160 CONTINUE 2060 END IF 2061 END IF 2062 ELSE 2063 IF( LSAME( TRANSA, 'N' ) )THEN 2064* 2065* Form B := alpha*B*inv( A ). 2066* 2067 IF( UPPER )THEN 2068 DO 210, J = 1, N 2069 IF( ALPHA.NE.ONE )THEN 2070 DO 170, I = 1, M 2071 B( I, J ) = ALPHA*B( I, J ) 2072 170 CONTINUE 2073 END IF 2074 DO 190, K = 1, J - 1 2075 IF( A( K, J ).NE.ZERO )THEN 2076 DO 180, I = 1, M 2077 B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 2078 180 CONTINUE 2079 END IF 2080 190 CONTINUE 2081 IF( NOUNIT )THEN 2082 TEMP = ONE/A( J, J ) 2083 DO 200, I = 1, M 2084 B( I, J ) = TEMP*B( I, J ) 2085 200 CONTINUE 2086 END IF 2087 210 CONTINUE 2088 ELSE 2089 DO 260, J = N, 1, -1 2090 IF( ALPHA.NE.ONE )THEN 2091 DO 220, I = 1, M 2092 B( I, J ) = ALPHA*B( I, J ) 2093 220 CONTINUE 2094 END IF 2095 DO 240, K = J + 1, N 2096 IF( A( K, J ).NE.ZERO )THEN 2097 DO 230, I = 1, M 2098 B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 2099 230 CONTINUE 2100 END IF 2101 240 CONTINUE 2102 IF( NOUNIT )THEN 2103 TEMP = ONE/A( J, J ) 2104 DO 250, I = 1, M 2105 B( I, J ) = TEMP*B( I, J ) 2106 250 CONTINUE 2107 END IF 2108 260 CONTINUE 2109 END IF 2110 ELSE 2111* 2112* Form B := alpha*B*inv( A' ). 2113* 2114 IF( UPPER )THEN 2115 DO 310, K = N, 1, -1 2116 IF( NOUNIT )THEN 2117 TEMP = ONE/A( K, K ) 2118 DO 270, I = 1, M 2119 B( I, K ) = TEMP*B( I, K ) 2120 270 CONTINUE 2121 END IF 2122 DO 290, J = 1, K - 1 2123 IF( A( J, K ).NE.ZERO )THEN 2124 TEMP = A( J, K ) 2125 DO 280, I = 1, M 2126 B( I, J ) = B( I, J ) - TEMP*B( I, K ) 2127 280 CONTINUE 2128 END IF 2129 290 CONTINUE 2130 IF( ALPHA.NE.ONE )THEN 2131 DO 300, I = 1, M 2132 B( I, K ) = ALPHA*B( I, K ) 2133 300 CONTINUE 2134 END IF 2135 310 CONTINUE 2136 ELSE 2137 DO 360, K = 1, N 2138 IF( NOUNIT )THEN 2139 TEMP = ONE/A( K, K ) 2140 DO 320, I = 1, M 2141 B( I, K ) = TEMP*B( I, K ) 2142 320 CONTINUE 2143 END IF 2144 DO 340, J = K + 1, N 2145 IF( A( J, K ).NE.ZERO )THEN 2146 TEMP = A( J, K ) 2147 DO 330, I = 1, M 2148 B( I, J ) = B( I, J ) - TEMP*B( I, K ) 2149 330 CONTINUE 2150 END IF 2151 340 CONTINUE 2152 IF( ALPHA.NE.ONE )THEN 2153 DO 350, I = 1, M 2154 B( I, K ) = ALPHA*B( I, K ) 2155 350 CONTINUE 2156 END IF 2157 360 CONTINUE 2158 END IF 2159 END IF 2160 END IF 2161* 2162 RETURN 2163* 2164* End of DTRSM . 2165* 2166 END 2167cc 2168 SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, 2169 $ BETA, C, LDC ) 2170* .. Scalar Arguments .. 2171 CHARACTER*1 TRANSA, TRANSB 2172 INTEGER M, N, K, LDA, LDB, LDC 2173 DOUBLE PRECISION ALPHA, BETA 2174* .. Array Arguments .. 2175 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) 2176* .. 2177* 2178* Purpose 2179* ======= 2180* 2181* DGEMM performs one of the matrix-matrix operations 2182* 2183* C := alpha*op( A )*op( B ) + beta*C, 2184* 2185* where op( X ) is one of 2186* 2187* op( X ) = X or op( X ) = X', 2188* 2189* alpha and beta are scalars, and A, B and C are matrices, with op( A ) 2190* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. 2191* 2192* Parameters 2193* ========== 2194* 2195* TRANSA - CHARACTER*1. 2196* On entry, TRANSA specifies the form of op( A ) to be used in 2197* the matrix multiplication as follows: 2198* 2199* TRANSA = 'N' or 'n', op( A ) = A. 2200* 2201* TRANSA = 'T' or 't', op( A ) = A'. 2202* 2203* TRANSA = 'C' or 'c', op( A ) = A'. 2204* 2205* Unchanged on exit. 2206* 2207* TRANSB - CHARACTER*1. 2208* On entry, TRANSB specifies the form of op( B ) to be used in 2209* the matrix multiplication as follows: 2210* 2211* TRANSB = 'N' or 'n', op( B ) = B. 2212* 2213* TRANSB = 'T' or 't', op( B ) = B'. 2214* 2215* TRANSB = 'C' or 'c', op( B ) = B'. 2216* 2217* Unchanged on exit. 2218* 2219* M - INTEGER. 2220* On entry, M specifies the number of rows of the matrix 2221* op( A ) and of the matrix C. M must be at least zero. 2222* Unchanged on exit. 2223* 2224* N - INTEGER. 2225* On entry, N specifies the number of columns of the matrix 2226* op( B ) and the number of columns of the matrix C. N must be 2227* at least zero. 2228* Unchanged on exit. 2229* 2230* K - INTEGER. 2231* On entry, K specifies the number of columns of the matrix 2232* op( A ) and the number of rows of the matrix op( B ). K must 2233* be at least zero. 2234* Unchanged on exit. 2235* 2236* ALPHA - DOUBLE PRECISION. 2237* On entry, ALPHA specifies the scalar alpha. 2238* Unchanged on exit. 2239* 2240* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is 2241* k when TRANSA = 'N' or 'n', and is m otherwise. 2242* Before entry with TRANSA = 'N' or 'n', the leading m by k 2243* part of the array A must contain the matrix A, otherwise 2244* the leading k by m part of the array A must contain the 2245* matrix A. 2246* Unchanged on exit. 2247* 2248* LDA - INTEGER. 2249* On entry, LDA specifies the first dimension of A as declared 2250* in the calling (sub) program. When TRANSA = 'N' or 'n' then 2251* LDA must be at least max( 1, m ), otherwise LDA must be at 2252* least max( 1, k ). 2253* Unchanged on exit. 2254* 2255* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is 2256* n when TRANSB = 'N' or 'n', and is k otherwise. 2257* Before entry with TRANSB = 'N' or 'n', the leading k by n 2258* part of the array B must contain the matrix B, otherwise 2259* the leading n by k part of the array B must contain the 2260* matrix B. 2261* Unchanged on exit. 2262* 2263* LDB - INTEGER. 2264* On entry, LDB specifies the first dimension of B as declared 2265* in the calling (sub) program. When TRANSB = 'N' or 'n' then 2266* LDB must be at least max( 1, k ), otherwise LDB must be at 2267* least max( 1, n ). 2268* Unchanged on exit. 2269* 2270* BETA - DOUBLE PRECISION. 2271* On entry, BETA specifies the scalar beta. When BETA is 2272* supplied as zero then C need not be set on input. 2273* Unchanged on exit. 2274* 2275* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). 2276* Before entry, the leading m by n part of the array C must 2277* contain the matrix C, except when beta is zero, in which 2278* case C need not be set on entry. 2279* On exit, the array C is overwritten by the m by n matrix 2280* ( alpha*op( A )*op( B ) + beta*C ). 2281* 2282* LDC - INTEGER. 2283* On entry, LDC specifies the first dimension of C as declared 2284* in the calling (sub) program. LDC must be at least 2285* max( 1, m ). 2286* Unchanged on exit. 2287* 2288* 2289* Level 3 Blas routine. 2290* 2291* -- Written on 8-February-1989. 2292* Jack Dongarra, Argonne National Laboratory. 2293* Iain Duff, AERE Harwell. 2294* Jeremy Du Croz, Numerical Algorithms Group Ltd. 2295* Sven Hammarling, Numerical Algorithms Group Ltd. 2296* 2297* 2298* .. External Functions .. 2299 LOGICAL LSAME 2300 EXTERNAL LSAME 2301* .. External Subroutines .. 2302 EXTERNAL XERBLA 2303* .. Intrinsic Functions .. 2304 INTRINSIC MAX 2305* .. Local Scalars .. 2306 LOGICAL NOTA, NOTB 2307 INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB 2308 DOUBLE PRECISION TEMP 2309* .. Parameters .. 2310 DOUBLE PRECISION ONE , ZERO 2311 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 2312* .. 2313* .. Executable Statements .. 2314* 2315* Set NOTA and NOTB as true if A and B respectively are not 2316* transposed and set NROWA, NCOLA and NROWB as the number of rows 2317* and columns of A and the number of rows of B respectively. 2318* 2319 NOTA = LSAME( TRANSA, 'N' ) 2320 NOTB = LSAME( TRANSB, 'N' ) 2321 IF( NOTA )THEN 2322 NROWA = M 2323 NCOLA = K 2324 ELSE 2325 NROWA = K 2326 NCOLA = M 2327 END IF 2328 IF( NOTB )THEN 2329 NROWB = K 2330 ELSE 2331 NROWB = N 2332 END IF 2333* 2334* Test the input parameters. 2335* 2336 INFO = 0 2337 IF( ( .NOT.NOTA ).AND. 2338 $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. 2339 $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN 2340 INFO = 1 2341 ELSE IF( ( .NOT.NOTB ).AND. 2342 $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. 2343 $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN 2344 INFO = 2 2345 ELSE IF( M .LT.0 )THEN 2346 INFO = 3 2347 ELSE IF( N .LT.0 )THEN 2348 INFO = 4 2349 ELSE IF( K .LT.0 )THEN 2350 INFO = 5 2351 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN 2352 INFO = 8 2353 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN 2354 INFO = 10 2355 ELSE IF( LDC.LT.MAX( 1, M ) )THEN 2356 INFO = 13 2357 END IF 2358 IF( INFO.NE.0 )THEN 2359 CALL XERBLA( 'DGEMM ', INFO ) 2360 RETURN 2361 END IF 2362* 2363* Quick return if possible. 2364* 2365 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. 2366 $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) 2367 $ RETURN 2368* 2369* And if alpha.eq.zero. 2370* 2371 IF( ALPHA.EQ.ZERO )THEN 2372 IF( BETA.EQ.ZERO )THEN 2373 DO 20, J = 1, N 2374 DO 10, I = 1, M 2375 C( I, J ) = ZERO 2376 10 CONTINUE 2377 20 CONTINUE 2378 ELSE 2379 DO 40, J = 1, N 2380 DO 30, I = 1, M 2381 C( I, J ) = BETA*C( I, J ) 2382 30 CONTINUE 2383 40 CONTINUE 2384 END IF 2385 RETURN 2386 END IF 2387* 2388* Start the operations. 2389* 2390 IF( NOTB )THEN 2391 IF( NOTA )THEN 2392* 2393* Form C := alpha*A*B + beta*C. 2394* 2395 DO 90, J = 1, N 2396 IF( BETA.EQ.ZERO )THEN 2397 DO 50, I = 1, M 2398 C( I, J ) = ZERO 2399 50 CONTINUE 2400 ELSE IF( BETA.NE.ONE )THEN 2401 DO 60, I = 1, M 2402 C( I, J ) = BETA*C( I, J ) 2403 60 CONTINUE 2404 END IF 2405 DO 80, L = 1, K 2406 IF( B( L, J ).NE.ZERO )THEN 2407 TEMP = ALPHA*B( L, J ) 2408 DO 70, I = 1, M 2409 C( I, J ) = C( I, J ) + TEMP*A( I, L ) 2410 70 CONTINUE 2411 END IF 2412 80 CONTINUE 2413 90 CONTINUE 2414 ELSE 2415* 2416* Form C := alpha*A'*B + beta*C 2417* 2418 DO 120, J = 1, N 2419 DO 110, I = 1, M 2420 TEMP = ZERO 2421 DO 100, L = 1, K 2422 TEMP = TEMP + A( L, I )*B( L, J ) 2423 100 CONTINUE 2424 IF( BETA.EQ.ZERO )THEN 2425 C( I, J ) = ALPHA*TEMP 2426 ELSE 2427 C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) 2428 END IF 2429 110 CONTINUE 2430 120 CONTINUE 2431 END IF 2432 ELSE 2433 IF( NOTA )THEN 2434* 2435* Form C := alpha*A*B' + beta*C 2436* 2437 DO 170, J = 1, N 2438 IF( BETA.EQ.ZERO )THEN 2439 DO 130, I = 1, M 2440 C( I, J ) = ZERO 2441 130 CONTINUE 2442 ELSE IF( BETA.NE.ONE )THEN 2443 DO 140, I = 1, M 2444 C( I, J ) = BETA*C( I, J ) 2445 140 CONTINUE 2446 END IF 2447 DO 160, L = 1, K 2448 IF( B( J, L ).NE.ZERO )THEN 2449 TEMP = ALPHA*B( J, L ) 2450 DO 150, I = 1, M 2451 C( I, J ) = C( I, J ) + TEMP*A( I, L ) 2452 150 CONTINUE 2453 END IF 2454 160 CONTINUE 2455 170 CONTINUE 2456 ELSE 2457* 2458* Form C := alpha*A'*B' + beta*C 2459* 2460 DO 200, J = 1, N 2461 DO 190, I = 1, M 2462 TEMP = ZERO 2463 DO 180, L = 1, K 2464 TEMP = TEMP + A( L, I )*B( J, L ) 2465 180 CONTINUE 2466 IF( BETA.EQ.ZERO )THEN 2467 C( I, J ) = ALPHA*TEMP 2468 ELSE 2469 C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) 2470 END IF 2471 190 CONTINUE 2472 200 CONTINUE 2473 END IF 2474 END IF 2475* 2476 RETURN 2477* 2478* End of DGEMM . 2479* 2480 END 2481