1 subroutine daxpy(n,da,dx,incx,dy,incy) 2c 3c constant times a vector plus a vector. 4c uses unrolled loops for increments equal to one. 5c jack dongarra, linpack, 3/11/78. 6c 7 double precision dx(1),dy(1),da 8 integer i,incx,incy,ix,iy,m,mp1,n 9c 10 if(n.le.0)return 11 if (da .eq. 0.0d0) return 12 if(incx.eq.1.and.incy.eq.1)go to 20 13c 14c code for unequal increments or equal increments 15c not equal to 1 16c 17 ix = 1 18 iy = 1 19 if(incx.lt.0)ix = (-n+1)*incx + 1 20 if(incy.lt.0)iy = (-n+1)*incy + 1 21 do 10 i = 1,n 22 dy(iy) = dy(iy) + da*dx(ix) 23 ix = ix + incx 24 iy = iy + incy 25 10 continue 26 return 27c 28c code for both increments equal to 1 29c 30c 31c clean-up loop 32c 33 20 m = mod(n,4) 34 if( m .eq. 0 ) go to 40 35 do 30 i = 1,m 36 dy(i) = dy(i) + da*dx(i) 37 30 continue 38 if( n .lt. 4 ) return 39 40 mp1 = m + 1 40 do 50 i = mp1,n,4 41 dy(i) = dy(i) + da*dx(i) 42 dy(i + 1) = dy(i + 1) + da*dx(i + 1) 43 dy(i + 2) = dy(i + 2) + da*dx(i + 2) 44 dy(i + 3) = dy(i + 3) + da*dx(i + 3) 45 50 continue 46 return 47 end 48 subroutine dcopy(n,dx,incx,dy,incy) 49c 50c copies a vector, x, to a vector, y. 51c uses unrolled loops for increments equal to one. 52c jack dongarra, linpack, 3/11/78. 53c 54 double precision dx(1),dy(1) 55 integer i,incx,incy,ix,iy,m,mp1,n 56c 57 if(n.le.0)return 58 if(incx.eq.1.and.incy.eq.1)go to 20 59c 60c code for unequal increments or equal increments 61c not equal to 1 62c 63 ix = 1 64 iy = 1 65 if(incx.lt.0)ix = (-n+1)*incx + 1 66 if(incy.lt.0)iy = (-n+1)*incy + 1 67 do 10 i = 1,n 68 dy(iy) = dx(ix) 69 ix = ix + incx 70 iy = iy + incy 71 10 continue 72 return 73c 74c code for both increments equal to 1 75c 76c 77c clean-up loop 78c 79 20 m = mod(n,7) 80 if( m .eq. 0 ) go to 40 81 do 30 i = 1,m 82 dy(i) = dx(i) 83 30 continue 84 if( n .lt. 7 ) return 85 40 mp1 = m + 1 86 do 50 i = mp1,n,7 87 dy(i) = dx(i) 88 dy(i + 1) = dx(i + 1) 89 dy(i + 2) = dx(i + 2) 90 dy(i + 3) = dx(i + 3) 91 dy(i + 4) = dx(i + 4) 92 dy(i + 5) = dx(i + 5) 93 dy(i + 6) = dx(i + 6) 94 50 continue 95 return 96 end 97 double precision function ddot(n,dx,incx,dy,incy) 98c 99c forms the dot product of two vectors. 100c uses unrolled loops for increments equal to one. 101c jack dongarra, linpack, 3/11/78. 102c 103 double precision dx(1),dy(1),dtemp 104 integer i,incx,incy,ix,iy,m,mp1,n 105c 106 ddot = 0.0d0 107 dtemp = 0.0d0 108 if(n.le.0)return 109 if(incx.eq.1.and.incy.eq.1)go to 20 110c 111c code for unequal increments or equal increments 112c not equal to 1 113c 114 ix = 1 115 iy = 1 116 if(incx.lt.0)ix = (-n+1)*incx + 1 117 if(incy.lt.0)iy = (-n+1)*incy + 1 118 do 10 i = 1,n 119 dtemp = dtemp + dx(ix)*dy(iy) 120 ix = ix + incx 121 iy = iy + incy 122 10 continue 123 ddot = dtemp 124 return 125c 126c code for both increments equal to 1 127c 128c 129c clean-up loop 130c 131 20 m = mod(n,5) 132 if( m .eq. 0 ) go to 40 133 do 30 i = 1,m 134 dtemp = dtemp + dx(i)*dy(i) 135 30 continue 136 if( n .lt. 5 ) go to 60 137 40 mp1 = m + 1 138 do 50 i = mp1,n,5 139 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + 140 * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 141 50 continue 142 60 ddot = dtemp 143 144 return 145 end 146 147 SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, 148 $ BETA, C, LDC ) 149* .. Scalar Arguments .. 150 CHARACTER*1 TRANSA, TRANSB 151 INTEGER M, N, K, LDA, LDB, LDC 152 DOUBLE PRECISION ALPHA, BETA 153* .. Array Arguments .. 154 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) 155* .. 156* 157* Purpose 158* ======= 159* 160* DGEMM performs one of the matrix-matrix operations 161* 162* C := alpha*op( A )*op( B ) + beta*C, 163* 164* where op( X ) is one of 165* 166* op( X ) = X or op( X ) = X', 167* 168* alpha and beta are scalars, and A, B and C are matrices, with op( A ) 169* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. 170* 171* Parameters 172* ========== 173* 174* TRANSA - CHARACTER*1. 175* On entry, TRANSA specifies the form of op( A ) to be used in 176* the matrix multiplication as follows: 177* 178* TRANSA = 'N' or 'n', op( A ) = A. 179* 180* TRANSA = 'T' or 't', op( A ) = A'. 181* 182* TRANSA = 'C' or 'c', op( A ) = A'. 183* 184* Unchanged on exit. 185* 186* TRANSB - CHARACTER*1. 187* On entry, TRANSB specifies the form of op( B ) to be used in 188* the matrix multiplication as follows: 189* 190* TRANSB = 'N' or 'n', op( B ) = B. 191* 192* TRANSB = 'T' or 't', op( B ) = B'. 193* 194* TRANSB = 'C' or 'c', op( B ) = B'. 195* 196* Unchanged on exit. 197* 198* M - INTEGER. 199* On entry, M specifies the number of rows of the matrix 200* op( A ) and of the matrix C. M must be at least zero. 201* Unchanged on exit. 202* 203* N - INTEGER. 204* On entry, N specifies the number of columns of the matrix 205* op( B ) and the number of columns of the matrix C. N must be 206* at least zero. 207* Unchanged on exit. 208* 209* K - INTEGER. 210* On entry, K specifies the number of columns of the matrix 211* op( A ) and the number of rows of the matrix op( B ). K must 212* be at least zero. 213* Unchanged on exit. 214* 215* ALPHA - DOUBLE PRECISION. 216* On entry, ALPHA specifies the scalar alpha. 217* Unchanged on exit. 218* 219* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is 220* k when TRANSA = 'N' or 'n', and is m otherwise. 221* Before entry with TRANSA = 'N' or 'n', the leading m by k 222* part of the array A must contain the matrix A, otherwise 223* the leading k by m part of the array A must contain the 224* matrix A. 225* Unchanged on exit. 226* 227* LDA - INTEGER. 228* On entry, LDA specifies the first dimension of A as declared 229* in the calling (sub) program. When TRANSA = 'N' or 'n' then 230* LDA must be at least max( 1, m ), otherwise LDA must be at 231* least max( 1, k ). 232* Unchanged on exit. 233* 234* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is 235* n when TRANSB = 'N' or 'n', and is k otherwise. 236* Before entry with TRANSB = 'N' or 'n', the leading k by n 237* part of the array B must contain the matrix B, otherwise 238* the leading n by k part of the array B must contain the 239* matrix B. 240* Unchanged on exit. 241* 242* LDB - INTEGER. 243* On entry, LDB specifies the first dimension of B as declared 244* in the calling (sub) program. When TRANSB = 'N' or 'n' then 245* LDB must be at least max( 1, k ), otherwise LDB must be at 246* least max( 1, n ). 247* Unchanged on exit. 248* 249* BETA - DOUBLE PRECISION. 250* On entry, BETA specifies the scalar beta. When BETA is 251* supplied as zero then C need not be set on input. 252* Unchanged on exit. 253* 254* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). 255* Before entry, the leading m by n part of the array C must 256* contain the matrix C, except when beta is zero, in which 257* case C need not be set on entry. 258* On exit, the array C is overwritten by the m by n matrix 259* ( alpha*op( A )*op( B ) + beta*C ). 260* 261* LDC - INTEGER. 262* On entry, LDC specifies the first dimension of C as declared 263* in the calling (sub) program. LDC must be at least 264* max( 1, m ). 265* Unchanged on exit. 266* 267* 268* Level 3 Blas routine. 269* 270* -- Written on 8-February-1989. 271* Jack Dongarra, Argonne National Laboratory. 272* Iain Duff, AERE Harwell. 273* Jeremy Du Croz, Numerical Algorithms Group Ltd. 274* Sven Hammarling, Numerical Algorithms Group Ltd. 275* 276* 277* .. External Functions .. 278 LOGICAL LSAME 279 EXTERNAL LSAME 280* .. External Subroutines .. 281 EXTERNAL XERBLA 282* .. Intrinsic Functions .. 283 INTRINSIC MAX 284* .. Local Scalars .. 285 LOGICAL NOTA, NOTB 286 INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB 287 DOUBLE PRECISION TEMP 288* .. Parameters .. 289 DOUBLE PRECISION ONE , ZERO 290 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 291* .. 292* .. Executable Statements .. 293* 294* Set NOTA and NOTB as true if A and B respectively are not 295* transposed and set NROWA, NCOLA and NROWB as the number of rows 296* and columns of A and the number of rows of B respectively. 297* 298 NOTA = LSAME( TRANSA, 'N' ) 299 NOTB = LSAME( TRANSB, 'N' ) 300 IF( NOTA )THEN 301 NROWA = M 302 NCOLA = K 303 ELSE 304 NROWA = K 305 NCOLA = M 306 END IF 307 IF( NOTB )THEN 308 NROWB = K 309 ELSE 310 NROWB = N 311 END IF 312* 313* Test the input parameters. 314* 315 INFO = 0 316 IF( ( .NOT.NOTA ).AND. 317 $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. 318 $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN 319 INFO = 1 320 ELSE IF( ( .NOT.NOTB ).AND. 321 $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. 322 $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN 323 INFO = 2 324 ELSE IF( M .LT.0 )THEN 325 INFO = 3 326 ELSE IF( N .LT.0 )THEN 327 INFO = 4 328 ELSE IF( K .LT.0 )THEN 329 INFO = 5 330 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN 331 INFO = 8 332 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN 333 INFO = 10 334 ELSE IF( LDC.LT.MAX( 1, M ) )THEN 335 INFO = 13 336 END IF 337 IF( INFO.NE.0 )THEN 338 CALL XERBLA( 'DGEMM ', INFO ) 339 RETURN 340 END IF 341* 342* Quick return if possible. 343* 344 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. 345 $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) 346 $ RETURN 347* 348* And if alpha.eq.zero. 349* 350 IF( ALPHA.EQ.ZERO )THEN 351 IF( BETA.EQ.ZERO )THEN 352 DO 20, J = 1, N 353 DO 10, I = 1, M 354 C( I, J ) = ZERO 355 10 CONTINUE 356 20 CONTINUE 357 ELSE 358 DO 40, J = 1, N 359 DO 30, I = 1, M 360 C( I, J ) = BETA*C( I, J ) 361 30 CONTINUE 362 40 CONTINUE 363 END IF 364 RETURN 365 END IF 366* 367* Start the operations. 368* 369 IF( NOTB )THEN 370 IF( NOTA )THEN 371* 372* Form C := alpha*A*B + beta*C. 373* 374 DO 90, J = 1, N 375 IF( BETA.EQ.ZERO )THEN 376 DO 50, I = 1, M 377 C( I, J ) = ZERO 378 50 CONTINUE 379 ELSE IF( BETA.NE.ONE )THEN 380 DO 60, I = 1, M 381 C( I, J ) = BETA*C( I, J ) 382 60 CONTINUE 383 END IF 384 DO 80, L = 1, K 385 IF( B( L, J ).NE.ZERO )THEN 386 TEMP = ALPHA*B( L, J ) 387 DO 70, I = 1, M 388 C( I, J ) = C( I, J ) + TEMP*A( I, L ) 389 70 CONTINUE 390 END IF 391 80 CONTINUE 392 90 CONTINUE 393 ELSE 394* 395* Form C := alpha*A'*B + beta*C 396* 397 DO 120, J = 1, N 398 DO 110, I = 1, M 399 TEMP = ZERO 400 DO 100, L = 1, K 401 TEMP = TEMP + A( L, I )*B( L, J ) 402 100 CONTINUE 403 IF( BETA.EQ.ZERO )THEN 404 C( I, J ) = ALPHA*TEMP 405 ELSE 406 C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) 407 END IF 408 110 CONTINUE 409 120 CONTINUE 410 END IF 411 ELSE 412 IF( NOTA )THEN 413* 414* Form C := alpha*A*B' + beta*C 415* 416 DO 170, J = 1, N 417 IF( BETA.EQ.ZERO )THEN 418 DO 130, I = 1, M 419 C( I, J ) = ZERO 420 130 CONTINUE 421 ELSE IF( BETA.NE.ONE )THEN 422 DO 140, I = 1, M 423 C( I, J ) = BETA*C( I, J ) 424 140 CONTINUE 425 END IF 426 DO 160, L = 1, K 427 IF( B( J, L ).NE.ZERO )THEN 428 TEMP = ALPHA*B( J, L ) 429 DO 150, I = 1, M 430 C( I, J ) = C( I, J ) + TEMP*A( I, L ) 431 150 CONTINUE 432 END IF 433 160 CONTINUE 434 170 CONTINUE 435 ELSE 436* 437* Form C := alpha*A'*B' + beta*C 438* 439 DO 200, J = 1, N 440 DO 190, I = 1, M 441 TEMP = ZERO 442 DO 180, L = 1, K 443 TEMP = TEMP + A( L, I )*B( J, L ) 444 180 CONTINUE 445 IF( BETA.EQ.ZERO )THEN 446 C( I, J ) = ALPHA*TEMP 447 ELSE 448 C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) 449 END IF 450 190 CONTINUE 451 200 CONTINUE 452 END IF 453 END IF 454* 455 RETURN 456* 457* End of DGEMM . 458* 459 END 460 SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, 461 $ BETA, Y, INCY ) 462* .. Scalar Arguments .. 463 DOUBLE PRECISION ALPHA, BETA 464 INTEGER INCX, INCY, LDA, M, N 465 CHARACTER*1 TRANS 466* .. Array Arguments .. 467 DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) 468* .. 469* 470* Purpose 471* ======= 472* 473* DGEMV performs one of the matrix-vector operations 474* 475* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, 476* 477* where alpha and beta are scalars, x and y are vectors and A is an 478* m by n matrix. 479* 480* Parameters 481* ========== 482* 483* TRANS - CHARACTER*1. 484* On entry, TRANS specifies the operation to be performed as 485* follows: 486* 487* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. 488* 489* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. 490* 491* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. 492* 493* Unchanged on exit. 494* 495* M - INTEGER. 496* On entry, M specifies the number of rows of the matrix A. 497* M must be at least zero. 498* Unchanged on exit. 499* 500* N - INTEGER. 501* On entry, N specifies the number of columns of the matrix A. 502* N must be at least zero. 503* Unchanged on exit. 504* 505* ALPHA - DOUBLE PRECISION. 506* On entry, ALPHA specifies the scalar alpha. 507* Unchanged on exit. 508* 509* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 510* Before entry, the leading m by n part of the array A must 511* contain the matrix of coefficients. 512* Unchanged on exit. 513* 514* LDA - INTEGER. 515* On entry, LDA specifies the first dimension of A as declared 516* in the calling (sub) program. LDA must be at least 517* max( 1, m ). 518* Unchanged on exit. 519* 520* X - DOUBLE PRECISION array of DIMENSION at least 521* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' 522* and at least 523* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. 524* Before entry, the incremented array X must contain the 525* vector x. 526* Unchanged on exit. 527* 528* INCX - INTEGER. 529* On entry, INCX specifies the increment for the elements of 530* X. INCX must not be zero. 531* Unchanged on exit. 532* 533* BETA - DOUBLE PRECISION. 534* On entry, BETA specifies the scalar beta. When BETA is 535* supplied as zero then Y need not be set on input. 536* Unchanged on exit. 537* 538* Y - DOUBLE PRECISION array of DIMENSION at least 539* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' 540* and at least 541* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. 542* Before entry with BETA non-zero, the incremented array Y 543* must contain the vector y. On exit, Y is overwritten by the 544* updated vector y. 545* 546* INCY - INTEGER. 547* On entry, INCY specifies the increment for the elements of 548* Y. INCY must not be zero. 549* Unchanged on exit. 550* 551* 552* Level 2 Blas routine. 553* 554* -- Written on 22-October-1986. 555* Jack Dongarra, Argonne National Lab. 556* Jeremy Du Croz, Nag Central Office. 557* Sven Hammarling, Nag Central Office. 558* Richard Hanson, Sandia National Labs. 559* 560* 561* .. Parameters .. 562 DOUBLE PRECISION ONE , ZERO 563 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 564* .. Local Scalars .. 565 DOUBLE PRECISION TEMP 566 INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY 567* .. External Functions .. 568 LOGICAL LSAME 569 EXTERNAL LSAME 570* .. External Subroutines .. 571 EXTERNAL XERBLA 572* .. Intrinsic Functions .. 573 INTRINSIC MAX 574* .. 575* .. Executable Statements .. 576* 577* Test the input parameters. 578* 579 INFO = 0 580 IF ( .NOT.LSAME( TRANS, 'N' ).AND. 581 $ .NOT.LSAME( TRANS, 'T' ).AND. 582 $ .NOT.LSAME( TRANS, 'C' ) )THEN 583 INFO = 1 584 ELSE IF( M.LT.0 )THEN 585 INFO = 2 586 ELSE IF( N.LT.0 )THEN 587 INFO = 3 588 ELSE IF( LDA.LT.MAX( 1, M ) )THEN 589 INFO = 6 590 ELSE IF( INCX.EQ.0 )THEN 591 INFO = 8 592 ELSE IF( INCY.EQ.0 )THEN 593 INFO = 11 594 END IF 595 IF( INFO.NE.0 )THEN 596 CALL XERBLA( 'DGEMV ', INFO ) 597 RETURN 598 END IF 599* 600* Quick return if possible. 601* 602 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. 603 $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) 604 $ RETURN 605* 606* Set LENX and LENY, the lengths of the vectors x and y, and set 607* up the start points in X and Y. 608* 609 IF( LSAME( TRANS, 'N' ) )THEN 610 LENX = N 611 LENY = M 612 ELSE 613 LENX = M 614 LENY = N 615 END IF 616 IF( INCX.GT.0 )THEN 617 KX = 1 618 ELSE 619 KX = 1 - ( LENX - 1 )*INCX 620 END IF 621 IF( INCY.GT.0 )THEN 622 KY = 1 623 ELSE 624 KY = 1 - ( LENY - 1 )*INCY 625 END IF 626* 627* Start the operations. In this version the elements of A are 628* accessed sequentially with one pass through A. 629* 630* First form y := beta*y. 631* 632 IF( BETA.NE.ONE )THEN 633 IF( INCY.EQ.1 )THEN 634 IF( BETA.EQ.ZERO )THEN 635 DO 10, I = 1, LENY 636 Y( I ) = ZERO 637 10 CONTINUE 638 ELSE 639 DO 20, I = 1, LENY 640 Y( I ) = BETA*Y( I ) 641 20 CONTINUE 642 END IF 643 ELSE 644 IY = KY 645 IF( BETA.EQ.ZERO )THEN 646 DO 30, I = 1, LENY 647 Y( IY ) = ZERO 648 IY = IY + INCY 649 30 CONTINUE 650 ELSE 651 DO 40, I = 1, LENY 652 Y( IY ) = BETA*Y( IY ) 653 IY = IY + INCY 654 40 CONTINUE 655 END IF 656 END IF 657 END IF 658 IF( ALPHA.EQ.ZERO ) 659 $ RETURN 660 IF( LSAME( TRANS, 'N' ) )THEN 661* 662* Form y := alpha*A*x + y. 663* 664 JX = KX 665 IF( INCY.EQ.1 )THEN 666 DO 60, J = 1, N 667 IF( X( JX ).NE.ZERO )THEN 668 TEMP = ALPHA*X( JX ) 669 DO 50, I = 1, M 670 Y( I ) = Y( I ) + TEMP*A( I, J ) 671 50 CONTINUE 672 END IF 673 JX = JX + INCX 674 60 CONTINUE 675 ELSE 676 DO 80, J = 1, N 677 IF( X( JX ).NE.ZERO )THEN 678 TEMP = ALPHA*X( JX ) 679 IY = KY 680 DO 70, I = 1, M 681 Y( IY ) = Y( IY ) + TEMP*A( I, J ) 682 IY = IY + INCY 683 70 CONTINUE 684 END IF 685 JX = JX + INCX 686 80 CONTINUE 687 END IF 688 ELSE 689* 690* Form y := alpha*A'*x + y. 691* 692 JY = KY 693 IF( INCX.EQ.1 )THEN 694 DO 100, J = 1, N 695 TEMP = ZERO 696 DO 90, I = 1, M 697 TEMP = TEMP + A( I, J )*X( I ) 698 90 CONTINUE 699 Y( JY ) = Y( JY ) + ALPHA*TEMP 700 JY = JY + INCY 701 100 CONTINUE 702 ELSE 703 DO 120, J = 1, N 704 TEMP = ZERO 705 IX = KX 706 DO 110, I = 1, M 707 TEMP = TEMP + A( I, J )*X( IX ) 708 IX = IX + INCX 709 110 CONTINUE 710 Y( JY ) = Y( JY ) + ALPHA*TEMP 711 JY = JY + INCY 712 120 CONTINUE 713 END IF 714 END IF 715* 716 RETURN 717* 718* End of DGEMV . 719* 720 END 721 SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) 722* .. Scalar Arguments .. 723 DOUBLE PRECISION ALPHA 724 INTEGER INCX, INCY, LDA, M, N 725* .. Array Arguments .. 726 DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) 727* .. 728* 729* Purpose 730* ======= 731* 732* DGER performs the rank 1 operation 733* 734* A := alpha*x*y' + A, 735* 736* where alpha is a scalar, x is an m element vector, y is an n element 737* vector and A is an m by n matrix. 738* 739* Parameters 740* ========== 741* 742* M - INTEGER. 743* On entry, M specifies the number of rows of the matrix A. 744* M must be at least zero. 745* Unchanged on exit. 746* 747* N - INTEGER. 748* On entry, N specifies the number of columns of the matrix A. 749* N must be at least zero. 750* Unchanged on exit. 751* 752* ALPHA - DOUBLE PRECISION. 753* On entry, ALPHA specifies the scalar alpha. 754* Unchanged on exit. 755* 756* X - DOUBLE PRECISION array of dimension at least 757* ( 1 + ( m - 1 )*abs( INCX ) ). 758* Before entry, the incremented array X must contain the m 759* element vector x. 760* Unchanged on exit. 761* 762* INCX - INTEGER. 763* On entry, INCX specifies the increment for the elements of 764* X. INCX must not be zero. 765* Unchanged on exit. 766* 767* Y - DOUBLE PRECISION array of dimension at least 768* ( 1 + ( n - 1 )*abs( INCY ) ). 769* Before entry, the incremented array Y must contain the n 770* element vector y. 771* Unchanged on exit. 772* 773* INCY - INTEGER. 774* On entry, INCY specifies the increment for the elements of 775* Y. INCY must not be zero. 776* Unchanged on exit. 777* 778* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 779* Before entry, the leading m by n part of the array A must 780* contain the matrix of coefficients. On exit, A is 781* overwritten by the updated matrix. 782* 783* LDA - INTEGER. 784* On entry, LDA specifies the first dimension of A as declared 785* in the calling (sub) program. LDA must be at least 786* max( 1, m ). 787* Unchanged on exit. 788* 789* 790* Level 2 Blas routine. 791* 792* -- Written on 22-October-1986. 793* Jack Dongarra, Argonne National Lab. 794* Jeremy Du Croz, Nag Central Office. 795* Sven Hammarling, Nag Central Office. 796* Richard Hanson, Sandia National Labs. 797* 798* 799* .. Parameters .. 800 DOUBLE PRECISION ZERO 801 PARAMETER ( ZERO = 0.0D+0 ) 802* .. Local Scalars .. 803 DOUBLE PRECISION TEMP 804 INTEGER I, INFO, IX, J, JY, KX 805* .. External Subroutines .. 806 EXTERNAL XERBLA 807* .. Intrinsic Functions .. 808 INTRINSIC MAX 809* .. 810* .. Executable Statements .. 811* 812* Test the input parameters. 813* 814 INFO = 0 815 IF ( M.LT.0 )THEN 816 INFO = 1 817 ELSE IF( N.LT.0 )THEN 818 INFO = 2 819 ELSE IF( INCX.EQ.0 )THEN 820 INFO = 5 821 ELSE IF( INCY.EQ.0 )THEN 822 INFO = 7 823 ELSE IF( LDA.LT.MAX( 1, M ) )THEN 824 INFO = 9 825 END IF 826 IF( INFO.NE.0 )THEN 827 CALL XERBLA( 'DGER ', INFO ) 828 RETURN 829 END IF 830* 831* Quick return if possible. 832* 833 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) 834 $ RETURN 835* 836* Start the operations. In this version the elements of A are 837* accessed sequentially with one pass through A. 838* 839 IF( INCY.GT.0 )THEN 840 JY = 1 841 ELSE 842 JY = 1 - ( N - 1 )*INCY 843 END IF 844 IF( INCX.EQ.1 )THEN 845 DO 20, J = 1, N 846 IF( Y( JY ).NE.ZERO )THEN 847 TEMP = ALPHA*Y( JY ) 848 DO 10, I = 1, M 849 A( I, J ) = A( I, J ) + X( I )*TEMP 850 10 CONTINUE 851 END IF 852 JY = JY + INCY 853 20 CONTINUE 854 ELSE 855 IF( INCX.GT.0 )THEN 856 KX = 1 857 ELSE 858 KX = 1 - ( M - 1 )*INCX 859 END IF 860 DO 40, J = 1, N 861 IF( Y( JY ).NE.ZERO )THEN 862 TEMP = ALPHA*Y( JY ) 863 IX = KX 864 DO 30, I = 1, M 865 A( I, J ) = A( I, J ) + X( IX )*TEMP 866 IX = IX + INCX 867 30 CONTINUE 868 END IF 869 JY = JY + INCY 870 40 CONTINUE 871 END IF 872* 873 RETURN 874* 875* End of DGER . 876* 877 END 878 double precision function dnrm2 ( n, dx, incx) 879 integer i, incx, ix, j, n, next 880 double precision dx(1), cutlo, cuthi, hitest, sum, xmax,zero,one 881 data zero, one /0.0d0, 1.0d0/ 882c 883c euclidean norm of the n-vector stored in dx() with storage 884c increment incx . 885c if n .le. 0 return with result = 0. 886c if n .ge. 1 then incx must be .ge. 1 887c 888c c.l.lawson, 1978 jan 08 889c modified to correct failure to update ix, 1/25/92. 890c modified 3/93 to return if incx .le. 0. 891c 892c four phase method using two built-in constants that are 893c hopefully applicable to all machines. 894c cutlo = maximum of dsqrt(u/eps) over all known machines. 895c cuthi = minimum of dsqrt(v) over all known machines. 896c where 897c eps = smallest no. such that eps + 1. .gt. 1. 898c u = smallest positive no. (underflow limit) 899c v = largest no. (overflow limit) 900c 901c brief outline of algorithm.. 902c 903c phase 1 scans zero components. 904c move to phase 2 when a component is nonzero and .le. cutlo 905c move to phase 3 when a component is .gt. cutlo 906c move to phase 4 when a component is .ge. cuthi/m 907c where m = n for x() real and m = 2*n for complex. 908c 909c values for cutlo and cuthi.. 910c from the environmental parameters listed in the imsl converter 911c document the limiting values are as follows.. 912c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are 913c univac and dec at 2**(-103) 914c thus cutlo = 2**(-51) = 4.44089e-16 915c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. 916c thus cuthi = 2**(63.5) = 1.30438e19 917c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. 918c thus cutlo = 2**(-33.5) = 8.23181d-11 919c cuthi, d.p. same as s.p. cuthi = 1.30438d19 920c data cutlo, cuthi / 8.232d-11, 1.304d19 / 921c data cutlo, cuthi / 4.441e-16, 1.304e19 / 922 data cutlo, cuthi / 8.232d-11, 1.304d19 / 923c 924 if(n .gt. 0 .and. incx.gt.0) go to 10 925 dnrm2 = zero 926 go to 300 927c 928 10 assign 30 to next 929 sum = zero 930 i = 1 931 ix = 1 932c begin main loop 933 20 go to next,(30, 50, 70, 110) 934 30 if( dabs(dx(i)) .gt. cutlo) go to 85 935 assign 50 to next 936 xmax = zero 937c 938c phase 1. sum is zero 939c 940 50 if( dx(i) .eq. zero) go to 200 941 if( dabs(dx(i)) .gt. cutlo) go to 85 942c 943c prepare for phase 2. 944 assign 70 to next 945 go to 105 946c 947c prepare for phase 4. 948c 949 100 continue 950 ix = j 951 assign 110 to next 952 sum = (sum / dx(i)) / dx(i) 953 105 xmax = dabs(dx(i)) 954 go to 115 955c 956c phase 2. sum is small. 957c scale to avoid destructive underflow. 958c 959 70 if( dabs(dx(i)) .gt. cutlo ) go to 75 960c 961c common code for phases 2 and 4. 962c in phase 4 sum is large. scale to avoid overflow. 963c 964 110 if( dabs(dx(i)) .le. xmax ) go to 115 965 sum = one + sum * (xmax / dx(i))**2 966 xmax = dabs(dx(i)) 967 go to 200 968c 969 115 sum = sum + (dx(i)/xmax)**2 970 go to 200 971c 972c 973c prepare for phase 3. 974c 975 75 sum = (sum * xmax) * xmax 976c 977c 978c for real or d.p. set hitest = cuthi/n 979c for complex set hitest = cuthi/(2*n) 980c 981 85 hitest = cuthi/float( n ) 982c 983c phase 3. sum is mid-range. no scaling. 984c 985 do 95 j = ix,n 986 if(dabs(dx(i)) .ge. hitest) go to 100 987 sum = sum + dx(i)**2 988 i = i + incx 989 95 continue 990 dnrm2 = dsqrt( sum ) 991 go to 300 992c 993 200 continue 994 ix = ix + 1 995 i = i + incx 996 if( ix .le. n ) go to 20 997c 998c end of main loop. 999c 1000c compute square root and adjust for scaling. 1001c 1002 dnrm2 = xmax * dsqrt(sum) 1003 300 continue 1004 1005 return 1006 end 1007 1008 subroutine drotg(da,db,c,s) 1009c 1010c construct givens plane rotation. 1011c jack dongarra, linpack, 3/11/78. 1012c 1013 double precision da,db,c,s,roe,scale,r,z 1014c 1015 roe = db 1016 if( dabs(da) .gt. dabs(db) ) roe = da 1017 scale = dabs(da) + dabs(db) 1018 if( scale .ne. 0.0d0 ) go to 10 1019 c = 1.0d0 1020 s = 0.0d0 1021 r = 0.0d0 1022 z = 0.0d0 1023 go to 20 1024 10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2) 1025 r = dsign(1.0d0,roe)*r 1026 c = da/r 1027 s = db/r 1028 z = 1.0d0 1029 if( dabs(da) .gt. dabs(db) ) z = s 1030 if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c 1031 20 da = r 1032 db = z 1033 return 1034 end 1035 subroutine drot (n,dx,incx,dy,incy,c,s) 1036c 1037c applies a plane rotation. 1038c jack dongarra, linpack, 3/11/78. 1039c 1040 double precision dx(1),dy(1),dtemp,c,s 1041 integer i,incx,incy,ix,iy,n 1042c 1043 if(n.le.0)return 1044 if(incx.eq.1.and.incy.eq.1)go to 20 1045c 1046c code for unequal increments or equal increments not equal 1047c to 1 1048c 1049 ix = 1 1050 iy = 1 1051 if(incx.lt.0)ix = (-n+1)*incx + 1 1052 if(incy.lt.0)iy = (-n+1)*incy + 1 1053 do 10 i = 1,n 1054 dtemp = c*dx(ix) + s*dy(iy) 1055 dy(iy) = c*dy(iy) - s*dx(ix) 1056 dx(ix) = dtemp 1057 ix = ix + incx 1058 iy = iy + incy 1059 10 continue 1060 return 1061c 1062c code for both increments equal to 1 1063c 1064 20 do 30 i = 1,n 1065 dtemp = c*dx(i) + s*dy(i) 1066 dy(i) = c*dy(i) - s*dx(i) 1067 dx(i) = dtemp 1068 30 continue 1069 return 1070 end 1071 subroutine dscal(n,da,dx,incx) 1072c 1073c scales a vector by a constant. 1074c uses unrolled loops for increment equal to one. 1075c jack dongarra, linpack, 3/11/78. 1076c modified 3/93 to return if incx .le. 0. 1077c 1078 double precision da,dx(1) 1079 integer i,incx,m,mp1,n,nincx 1080c 1081 if( n.le.0 .or. incx.le.0 )return 1082 if(incx.eq.1)go to 20 1083c 1084c code for increment not equal to 1 1085c 1086 nincx = n*incx 1087 do 10 i = 1,nincx,incx 1088 dx(i) = da*dx(i) 1089 10 continue 1090 return 1091c 1092c code for increment equal to 1 1093c 1094c 1095c clean-up loop 1096c 1097 20 m = mod(n,5) 1098 if( m .eq. 0 ) go to 40 1099 do 30 i = 1,m 1100 dx(i) = da*dx(i) 1101 30 continue 1102 if( n .lt. 5 ) return 1103 40 mp1 = m + 1 1104 do 50 i = mp1,n,5 1105 dx(i) = da*dx(i) 1106 dx(i + 1) = da*dx(i + 1) 1107 dx(i + 2) = da*dx(i + 2) 1108 dx(i + 3) = da*dx(i + 3) 1109 dx(i + 4) = da*dx(i + 4) 1110 50 continue 1111 return 1112 end 1113 subroutine dswap (n,dx,incx,dy,incy) 1114c 1115c interchanges two vectors. 1116c uses unrolled loops for increments equal one. 1117c jack dongarra, linpack, 3/11/78. 1118c 1119 double precision dx(1),dy(1),dtemp 1120 integer i,incx,incy,ix,iy,m,mp1,n 1121c 1122 if(n.le.0)return 1123 if(incx.eq.1.and.incy.eq.1)go to 20 1124c 1125c code for unequal increments or equal increments not equal 1126c to 1 1127c 1128 ix = 1 1129 iy = 1 1130 if(incx.lt.0)ix = (-n+1)*incx + 1 1131 if(incy.lt.0)iy = (-n+1)*incy + 1 1132 do 10 i = 1,n 1133 dtemp = dx(ix) 1134 dx(ix) = dy(iy) 1135 dy(iy) = dtemp 1136 ix = ix + incx 1137 iy = iy + incy 1138 10 continue 1139 return 1140c 1141c code for both increments equal to 1 1142c 1143c 1144c clean-up loop 1145c 1146 20 m = mod(n,3) 1147 if( m .eq. 0 ) go to 40 1148 do 30 i = 1,m 1149 dtemp = dx(i) 1150 dx(i) = dy(i) 1151 dy(i) = dtemp 1152 30 continue 1153 if( n .lt. 3 ) return 1154 40 mp1 = m + 1 1155 do 50 i = mp1,n,3 1156 dtemp = dx(i) 1157 dx(i) = dy(i) 1158 dy(i) = dtemp 1159 dtemp = dx(i + 1) 1160 dx(i + 1) = dy(i + 1) 1161 dy(i + 1) = dtemp 1162 dtemp = dx(i + 2) 1163 dx(i + 2) = dy(i + 2) 1164 dy(i + 2) = dtemp 1165 50 continue 1166 return 1167 end 1168 SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX, 1169 $ BETA, Y, INCY ) 1170* .. Scalar Arguments .. 1171 DOUBLE PRECISION ALPHA, BETA 1172 INTEGER INCX, INCY, LDA, N 1173 CHARACTER*1 UPLO 1174* .. Array Arguments .. 1175 DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) 1176* .. 1177* 1178* Purpose 1179* ======= 1180* 1181* DSYMV performs the matrix-vector operation 1182* 1183* y := alpha*A*x + beta*y, 1184* 1185* where alpha and beta are scalars, x and y are n element vectors and 1186* A is an n by n symmetric matrix. 1187* 1188* Parameters 1189* ========== 1190* 1191* UPLO - CHARACTER*1. 1192* On entry, UPLO specifies whether the upper or lower 1193* triangular part of the array A is to be referenced as 1194* follows: 1195* 1196* UPLO = 'U' or 'u' Only the upper triangular part of A 1197* is to be referenced. 1198* 1199* UPLO = 'L' or 'l' Only the lower triangular part of A 1200* is to be referenced. 1201* 1202* Unchanged on exit. 1203* 1204* N - INTEGER. 1205* On entry, N specifies the order of the matrix A. 1206* N must be at least zero. 1207* Unchanged on exit. 1208* 1209* ALPHA - DOUBLE PRECISION. 1210* On entry, ALPHA specifies the scalar alpha. 1211* Unchanged on exit. 1212* 1213* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 1214* Before entry with UPLO = 'U' or 'u', the leading n by n 1215* upper triangular part of the array A must contain the upper 1216* triangular part of the symmetric matrix and the strictly 1217* lower triangular part of A is not referenced. 1218* Before entry with UPLO = 'L' or 'l', the leading n by n 1219* lower triangular part of the array A must contain the lower 1220* triangular part of the symmetric matrix and the strictly 1221* upper triangular part of A is not referenced. 1222* Unchanged on exit. 1223* 1224* LDA - INTEGER. 1225* On entry, LDA specifies the first dimension of A as declared 1226* in the calling (sub) program. LDA must be at least 1227* max( 1, n ). 1228* Unchanged on exit. 1229* 1230* X - DOUBLE PRECISION array of dimension at least 1231* ( 1 + ( n - 1 )*abs( INCX ) ). 1232* Before entry, the incremented array X must contain the n 1233* element vector x. 1234* Unchanged on exit. 1235* 1236* INCX - INTEGER. 1237* On entry, INCX specifies the increment for the elements of 1238* X. INCX must not be zero. 1239* Unchanged on exit. 1240* 1241* BETA - DOUBLE PRECISION. 1242* On entry, BETA specifies the scalar beta. When BETA is 1243* supplied as zero then Y need not be set on input. 1244* Unchanged on exit. 1245* 1246* Y - DOUBLE PRECISION array of dimension at least 1247* ( 1 + ( n - 1 )*abs( INCY ) ). 1248* Before entry, the incremented array Y must contain the n 1249* element vector y. On exit, Y is overwritten by the updated 1250* vector y. 1251* 1252* INCY - INTEGER. 1253* On entry, INCY specifies the increment for the elements of 1254* Y. INCY must not be zero. 1255* Unchanged on exit. 1256* 1257* 1258* Level 2 Blas routine. 1259* 1260* -- Written on 22-October-1986. 1261* Jack Dongarra, Argonne National Lab. 1262* Jeremy Du Croz, Nag Central Office. 1263* Sven Hammarling, Nag Central Office. 1264* Richard Hanson, Sandia National Labs. 1265* 1266* 1267* .. Parameters .. 1268 DOUBLE PRECISION ONE , ZERO 1269 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 1270* .. Local Scalars .. 1271 DOUBLE PRECISION TEMP1, TEMP2 1272 INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY 1273* .. External Functions .. 1274 LOGICAL LSAME 1275 EXTERNAL LSAME 1276* .. External Subroutines .. 1277 EXTERNAL XERBLA 1278* .. Intrinsic Functions .. 1279 INTRINSIC MAX 1280* .. 1281* .. Executable Statements .. 1282* 1283* Test the input parameters. 1284* 1285 INFO = 0 1286 IF ( .NOT.LSAME( UPLO, 'U' ).AND. 1287 $ .NOT.LSAME( UPLO, 'L' ) )THEN 1288 INFO = 1 1289 ELSE IF( N.LT.0 )THEN 1290 INFO = 2 1291 ELSE IF( LDA.LT.MAX( 1, N ) )THEN 1292 INFO = 5 1293 ELSE IF( INCX.EQ.0 )THEN 1294 INFO = 7 1295 ELSE IF( INCY.EQ.0 )THEN 1296 INFO = 10 1297 END IF 1298 IF( INFO.NE.0 )THEN 1299 CALL XERBLA( 'DSYMV ', INFO ) 1300 RETURN 1301 END IF 1302* 1303* Quick return if possible. 1304* 1305 IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) 1306 $ RETURN 1307* 1308* Set up the start points in X and Y. 1309* 1310 IF( INCX.GT.0 )THEN 1311 KX = 1 1312 ELSE 1313 KX = 1 - ( N - 1 )*INCX 1314 END IF 1315 IF( INCY.GT.0 )THEN 1316 KY = 1 1317 ELSE 1318 KY = 1 - ( N - 1 )*INCY 1319 END IF 1320* 1321* Start the operations. In this version the elements of A are 1322* accessed sequentially with one pass through the triangular part 1323* of A. 1324* 1325* First form y := beta*y. 1326* 1327 IF( BETA.NE.ONE )THEN 1328 IF( INCY.EQ.1 )THEN 1329 IF( BETA.EQ.ZERO )THEN 1330 DO 10, I = 1, N 1331 Y( I ) = ZERO 1332 10 CONTINUE 1333 ELSE 1334 DO 20, I = 1, N 1335 Y( I ) = BETA*Y( I ) 1336 20 CONTINUE 1337 END IF 1338 ELSE 1339 IY = KY 1340 IF( BETA.EQ.ZERO )THEN 1341 DO 30, I = 1, N 1342 Y( IY ) = ZERO 1343 IY = IY + INCY 1344 30 CONTINUE 1345 ELSE 1346 DO 40, I = 1, N 1347 Y( IY ) = BETA*Y( IY ) 1348 IY = IY + INCY 1349 40 CONTINUE 1350 END IF 1351 END IF 1352 END IF 1353 IF( ALPHA.EQ.ZERO ) 1354 $ RETURN 1355 IF( LSAME( UPLO, 'U' ) )THEN 1356* 1357* Form y when A is stored in upper triangle. 1358* 1359 IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN 1360 DO 60, J = 1, N 1361 TEMP1 = ALPHA*X( J ) 1362 TEMP2 = ZERO 1363 DO 50, I = 1, J - 1 1364 Y( I ) = Y( I ) + TEMP1*A( I, J ) 1365 TEMP2 = TEMP2 + A( I, J )*X( I ) 1366 50 CONTINUE 1367 Y( J ) = Y( J ) + TEMP1*A( J, J ) + ALPHA*TEMP2 1368 60 CONTINUE 1369 ELSE 1370 JX = KX 1371 JY = KY 1372 DO 80, J = 1, N 1373 TEMP1 = ALPHA*X( JX ) 1374 TEMP2 = ZERO 1375 IX = KX 1376 IY = KY 1377 DO 70, I = 1, J - 1 1378 Y( IY ) = Y( IY ) + TEMP1*A( I, J ) 1379 TEMP2 = TEMP2 + A( I, J )*X( IX ) 1380 IX = IX + INCX 1381 IY = IY + INCY 1382 70 CONTINUE 1383 Y( JY ) = Y( JY ) + TEMP1*A( J, J ) + ALPHA*TEMP2 1384 JX = JX + INCX 1385 JY = JY + INCY 1386 80 CONTINUE 1387 END IF 1388 ELSE 1389* 1390* Form y when A is stored in lower triangle. 1391* 1392 IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN 1393 DO 100, J = 1, N 1394 TEMP1 = ALPHA*X( J ) 1395 TEMP2 = ZERO 1396 Y( J ) = Y( J ) + TEMP1*A( J, J ) 1397 DO 90, I = J + 1, N 1398 Y( I ) = Y( I ) + TEMP1*A( I, J ) 1399 TEMP2 = TEMP2 + A( I, J )*X( I ) 1400 90 CONTINUE 1401 Y( J ) = Y( J ) + ALPHA*TEMP2 1402 100 CONTINUE 1403 ELSE 1404 JX = KX 1405 JY = KY 1406 DO 120, J = 1, N 1407 TEMP1 = ALPHA*X( JX ) 1408 TEMP2 = ZERO 1409 Y( JY ) = Y( JY ) + TEMP1*A( J, J ) 1410 IX = JX 1411 IY = JY 1412 DO 110, I = J + 1, N 1413 IX = IX + INCX 1414 IY = IY + INCY 1415 Y( IY ) = Y( IY ) + TEMP1*A( I, J ) 1416 TEMP2 = TEMP2 + A( I, J )*X( IX ) 1417 110 CONTINUE 1418 Y( JY ) = Y( JY ) + ALPHA*TEMP2 1419 JX = JX + INCX 1420 JY = JY + INCY 1421 120 CONTINUE 1422 END IF 1423 END IF 1424* 1425 RETURN 1426* 1427* End of DSYMV . 1428* 1429 END 1430 SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) 1431* .. Scalar Arguments .. 1432 DOUBLE PRECISION ALPHA 1433 INTEGER INCX, INCY, LDA, N 1434 CHARACTER*1 UPLO 1435* .. Array Arguments .. 1436 DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) 1437* .. 1438* 1439* Purpose 1440* ======= 1441* 1442* DSYR2 performs the symmetric rank 2 operation 1443* 1444* A := alpha*x*y' + alpha*y*x' + A, 1445* 1446* where alpha is a scalar, x and y are n element vectors and A is an n 1447* by n symmetric matrix. 1448* 1449* Parameters 1450* ========== 1451* 1452* UPLO - CHARACTER*1. 1453* On entry, UPLO specifies whether the upper or lower 1454* triangular part of the array A is to be referenced as 1455* follows: 1456* 1457* UPLO = 'U' or 'u' Only the upper triangular part of A 1458* is to be referenced. 1459* 1460* UPLO = 'L' or 'l' Only the lower triangular part of A 1461* is to be referenced. 1462* 1463* Unchanged on exit. 1464* 1465* N - INTEGER. 1466* On entry, N specifies the order of the matrix A. 1467* N must be at least zero. 1468* Unchanged on exit. 1469* 1470* ALPHA - DOUBLE PRECISION. 1471* On entry, ALPHA specifies the scalar alpha. 1472* Unchanged on exit. 1473* 1474* X - DOUBLE PRECISION array of dimension at least 1475* ( 1 + ( n - 1 )*abs( INCX ) ). 1476* Before entry, the incremented array X must contain the n 1477* element vector x. 1478* Unchanged on exit. 1479* 1480* INCX - INTEGER. 1481* On entry, INCX specifies the increment for the elements of 1482* X. INCX must not be zero. 1483* Unchanged on exit. 1484* 1485* Y - DOUBLE PRECISION array of dimension at least 1486* ( 1 + ( n - 1 )*abs( INCY ) ). 1487* Before entry, the incremented array Y must contain the n 1488* element vector y. 1489* Unchanged on exit. 1490* 1491* INCY - INTEGER. 1492* On entry, INCY specifies the increment for the elements of 1493* Y. INCY must not be zero. 1494* Unchanged on exit. 1495* 1496* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 1497* Before entry with UPLO = 'U' or 'u', the leading n by n 1498* upper triangular part of the array A must contain the upper 1499* triangular part of the symmetric matrix and the strictly 1500* lower triangular part of A is not referenced. On exit, the 1501* upper triangular part of the array A is overwritten by the 1502* upper triangular part of the updated matrix. 1503* Before entry with UPLO = 'L' or 'l', the leading n by n 1504* lower triangular part of the array A must contain the lower 1505* triangular part of the symmetric matrix and the strictly 1506* upper triangular part of A is not referenced. On exit, the 1507* lower triangular part of the array A is overwritten by the 1508* lower triangular part of the updated matrix. 1509* 1510* LDA - INTEGER. 1511* On entry, LDA specifies the first dimension of A as declared 1512* in the calling (sub) program. LDA must be at least 1513* max( 1, n ). 1514* Unchanged on exit. 1515* 1516* 1517* Level 2 Blas routine. 1518* 1519* -- Written on 22-October-1986. 1520* Jack Dongarra, Argonne National Lab. 1521* Jeremy Du Croz, Nag Central Office. 1522* Sven Hammarling, Nag Central Office. 1523* Richard Hanson, Sandia National Labs. 1524* 1525* 1526* .. Parameters .. 1527 DOUBLE PRECISION ZERO 1528 PARAMETER ( ZERO = 0.0D+0 ) 1529* .. Local Scalars .. 1530 DOUBLE PRECISION TEMP1, TEMP2 1531 INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY 1532* .. External Functions .. 1533 LOGICAL LSAME 1534 EXTERNAL LSAME 1535* .. External Subroutines .. 1536 EXTERNAL XERBLA 1537* .. Intrinsic Functions .. 1538 INTRINSIC MAX 1539* .. 1540* .. Executable Statements .. 1541* 1542* Test the input parameters. 1543* 1544 INFO = 0 1545 IF ( .NOT.LSAME( UPLO, 'U' ).AND. 1546 $ .NOT.LSAME( UPLO, 'L' ) )THEN 1547 INFO = 1 1548 ELSE IF( N.LT.0 )THEN 1549 INFO = 2 1550 ELSE IF( INCX.EQ.0 )THEN 1551 INFO = 5 1552 ELSE IF( INCY.EQ.0 )THEN 1553 INFO = 7 1554 ELSE IF( LDA.LT.MAX( 1, N ) )THEN 1555 INFO = 9 1556 END IF 1557 IF( INFO.NE.0 )THEN 1558 CALL XERBLA( 'DSYR2 ', INFO ) 1559 RETURN 1560 END IF 1561* 1562* Quick return if possible. 1563* 1564 IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) 1565 $ RETURN 1566* 1567* Set up the start points in X and Y if the increments are not both 1568* unity. 1569* 1570 IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN 1571 IF( INCX.GT.0 )THEN 1572 KX = 1 1573 ELSE 1574 KX = 1 - ( N - 1 )*INCX 1575 END IF 1576 IF( INCY.GT.0 )THEN 1577 KY = 1 1578 ELSE 1579 KY = 1 - ( N - 1 )*INCY 1580 END IF 1581 JX = KX 1582 JY = KY 1583 END IF 1584* 1585* Start the operations. In this version the elements of A are 1586* accessed sequentially with one pass through the triangular part 1587* of A. 1588* 1589 IF( LSAME( UPLO, 'U' ) )THEN 1590* 1591* Form A when A is stored in the upper triangle. 1592* 1593 IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN 1594 DO 20, J = 1, N 1595 IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN 1596 TEMP1 = ALPHA*Y( J ) 1597 TEMP2 = ALPHA*X( J ) 1598 DO 10, I = 1, J 1599 A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 1600 10 CONTINUE 1601 END IF 1602 20 CONTINUE 1603 ELSE 1604 DO 40, J = 1, N 1605 IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN 1606 TEMP1 = ALPHA*Y( JY ) 1607 TEMP2 = ALPHA*X( JX ) 1608 IX = KX 1609 IY = KY 1610 DO 30, I = 1, J 1611 A( I, J ) = A( I, J ) + X( IX )*TEMP1 1612 $ + Y( IY )*TEMP2 1613 IX = IX + INCX 1614 IY = IY + INCY 1615 30 CONTINUE 1616 END IF 1617 JX = JX + INCX 1618 JY = JY + INCY 1619 40 CONTINUE 1620 END IF 1621 ELSE 1622* 1623* Form A when A is stored in the lower triangle. 1624* 1625 IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN 1626 DO 60, J = 1, N 1627 IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN 1628 TEMP1 = ALPHA*Y( J ) 1629 TEMP2 = ALPHA*X( J ) 1630 DO 50, I = J, N 1631 A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 1632 50 CONTINUE 1633 END IF 1634 60 CONTINUE 1635 ELSE 1636 DO 80, J = 1, N 1637 IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN 1638 TEMP1 = ALPHA*Y( JY ) 1639 TEMP2 = ALPHA*X( JX ) 1640 IX = JX 1641 IY = JY 1642 DO 70, I = J, N 1643 A( I, J ) = A( I, J ) + X( IX )*TEMP1 1644 $ + Y( IY )*TEMP2 1645 IX = IX + INCX 1646 IY = IY + INCY 1647 70 CONTINUE 1648 END IF 1649 JX = JX + INCX 1650 JY = JY + INCY 1651 80 CONTINUE 1652 END IF 1653 END IF 1654* 1655 RETURN 1656* 1657* End of DSYR2 . 1658* 1659 END 1660 SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, 1661 $ BETA, C, LDC ) 1662* .. Scalar Arguments .. 1663 CHARACTER*1 UPLO, TRANS 1664 INTEGER N, K, LDA, LDB, LDC 1665 DOUBLE PRECISION ALPHA, BETA 1666* .. Array Arguments .. 1667 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) 1668* .. 1669* 1670* Purpose 1671* ======= 1672* 1673* DSYR2K performs one of the symmetric rank 2k operations 1674* 1675* C := alpha*A*B' + alpha*B*A' + beta*C, 1676* 1677* or 1678* 1679* C := alpha*A'*B + alpha*B'*A + beta*C, 1680* 1681* where alpha and beta are scalars, C is an n by n symmetric matrix 1682* and A and B are n by k matrices in the first case and k by n 1683* matrices in the second case. 1684* 1685* Parameters 1686* ========== 1687* 1688* UPLO - CHARACTER*1. 1689* On entry, UPLO specifies whether the upper or lower 1690* triangular part of the array C is to be referenced as 1691* follows: 1692* 1693* UPLO = 'U' or 'u' Only the upper triangular part of C 1694* is to be referenced. 1695* 1696* UPLO = 'L' or 'l' Only the lower triangular part of C 1697* is to be referenced. 1698* 1699* Unchanged on exit. 1700* 1701* TRANS - CHARACTER*1. 1702* On entry, TRANS specifies the operation to be performed as 1703* follows: 1704* 1705* TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' + 1706* beta*C. 1707* 1708* TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A + 1709* beta*C. 1710* 1711* TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A + 1712* beta*C. 1713* 1714* Unchanged on exit. 1715* 1716* N - INTEGER. 1717* On entry, N specifies the order of the matrix C. N must be 1718* at least zero. 1719* Unchanged on exit. 1720* 1721* K - INTEGER. 1722* On entry with TRANS = 'N' or 'n', K specifies the number 1723* of columns of the matrices A and B, and on entry with 1724* TRANS = 'T' or 't' or 'C' or 'c', K specifies the number 1725* of rows of the matrices A and B. K must be at least zero. 1726* Unchanged on exit. 1727* 1728* ALPHA - DOUBLE PRECISION. 1729* On entry, ALPHA specifies the scalar alpha. 1730* Unchanged on exit. 1731* 1732* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is 1733* k when TRANS = 'N' or 'n', and is n otherwise. 1734* Before entry with TRANS = 'N' or 'n', the leading n by k 1735* part of the array A must contain the matrix A, otherwise 1736* the leading k by n part of the array A must contain the 1737* matrix A. 1738* Unchanged on exit. 1739* 1740* LDA - INTEGER. 1741* On entry, LDA specifies the first dimension of A as declared 1742* in the calling (sub) program. When TRANS = 'N' or 'n' 1743* then LDA must be at least max( 1, n ), otherwise LDA must 1744* be at least max( 1, k ). 1745* Unchanged on exit. 1746* 1747* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is 1748* k when TRANS = 'N' or 'n', and is n otherwise. 1749* Before entry with TRANS = 'N' or 'n', the leading n by k 1750* part of the array B must contain the matrix B, otherwise 1751* the leading k by n part of the array B must contain the 1752* matrix B. 1753* Unchanged on exit. 1754* 1755* LDB - INTEGER. 1756* On entry, LDB specifies the first dimension of B as declared 1757* in the calling (sub) program. When TRANS = 'N' or 'n' 1758* then LDB must be at least max( 1, n ), otherwise LDB must 1759* be at least max( 1, k ). 1760* Unchanged on exit. 1761* 1762* BETA - DOUBLE PRECISION. 1763* On entry, BETA specifies the scalar beta. 1764* Unchanged on exit. 1765* 1766* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). 1767* Before entry with UPLO = 'U' or 'u', the leading n by n 1768* upper triangular part of the array C must contain the upper 1769* triangular part of the symmetric matrix and the strictly 1770* lower triangular part of C is not referenced. On exit, the 1771* upper triangular part of the array C is overwritten by the 1772* upper triangular part of the updated matrix. 1773* Before entry with UPLO = 'L' or 'l', the leading n by n 1774* lower triangular part of the array C must contain the lower 1775* triangular part of the symmetric matrix and the strictly 1776* upper triangular part of C is not referenced. On exit, the 1777* lower triangular part of the array C is overwritten by the 1778* lower triangular part of the updated matrix. 1779* 1780* LDC - INTEGER. 1781* On entry, LDC specifies the first dimension of C as declared 1782* in the calling (sub) program. LDC must be at least 1783* max( 1, n ). 1784* Unchanged on exit. 1785* 1786* 1787* Level 3 Blas routine. 1788* 1789* 1790* -- Written on 8-February-1989. 1791* Jack Dongarra, Argonne National Laboratory. 1792* Iain Duff, AERE Harwell. 1793* Jeremy Du Croz, Numerical Algorithms Group Ltd. 1794* Sven Hammarling, Numerical Algorithms Group Ltd. 1795* 1796* 1797* .. External Functions .. 1798 LOGICAL LSAME 1799 EXTERNAL LSAME 1800* .. External Subroutines .. 1801 EXTERNAL XERBLA 1802* .. Intrinsic Functions .. 1803 INTRINSIC MAX 1804* .. Local Scalars .. 1805 LOGICAL UPPER 1806 INTEGER I, INFO, J, L, NROWA 1807 DOUBLE PRECISION TEMP1, TEMP2 1808* .. Parameters .. 1809 DOUBLE PRECISION ONE , ZERO 1810 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 1811* .. 1812* .. Executable Statements .. 1813* 1814* Test the input parameters. 1815* 1816 IF( LSAME( TRANS, 'N' ) )THEN 1817 NROWA = N 1818 ELSE 1819 NROWA = K 1820 END IF 1821 UPPER = LSAME( UPLO, 'U' ) 1822* 1823 INFO = 0 1824 IF( ( .NOT.UPPER ).AND. 1825 $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN 1826 INFO = 1 1827 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. 1828 $ ( .NOT.LSAME( TRANS, 'T' ) ).AND. 1829 $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN 1830 INFO = 2 1831 ELSE IF( N .LT.0 )THEN 1832 INFO = 3 1833 ELSE IF( K .LT.0 )THEN 1834 INFO = 4 1835 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN 1836 INFO = 7 1837 ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN 1838 INFO = 9 1839 ELSE IF( LDC.LT.MAX( 1, N ) )THEN 1840 INFO = 12 1841 END IF 1842 IF( INFO.NE.0 )THEN 1843 CALL XERBLA( 'DSYR2K', INFO ) 1844 RETURN 1845 END IF 1846* 1847* Quick return if possible. 1848* 1849 IF( ( N.EQ.0 ).OR. 1850 $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) 1851 $ RETURN 1852* 1853* And when alpha.eq.zero. 1854* 1855 IF( ALPHA.EQ.ZERO )THEN 1856 IF( UPPER )THEN 1857 IF( BETA.EQ.ZERO )THEN 1858 DO 20, J = 1, N 1859 DO 10, I = 1, J 1860 C( I, J ) = ZERO 1861 10 CONTINUE 1862 20 CONTINUE 1863 ELSE 1864 DO 40, J = 1, N 1865 DO 30, I = 1, J 1866 C( I, J ) = BETA*C( I, J ) 1867 30 CONTINUE 1868 40 CONTINUE 1869 END IF 1870 ELSE 1871 IF( BETA.EQ.ZERO )THEN 1872 DO 60, J = 1, N 1873 DO 50, I = J, N 1874 C( I, J ) = ZERO 1875 50 CONTINUE 1876 60 CONTINUE 1877 ELSE 1878 DO 80, J = 1, N 1879 DO 70, I = J, N 1880 C( I, J ) = BETA*C( I, J ) 1881 70 CONTINUE 1882 80 CONTINUE 1883 END IF 1884 END IF 1885 RETURN 1886 END IF 1887* 1888* Start the operations. 1889* 1890 IF( LSAME( TRANS, 'N' ) )THEN 1891* 1892* Form C := alpha*A*B' + alpha*B*A' + C. 1893* 1894 IF( UPPER )THEN 1895 DO 130, J = 1, N 1896 IF( BETA.EQ.ZERO )THEN 1897 DO 90, I = 1, J 1898 C( I, J ) = ZERO 1899 90 CONTINUE 1900 ELSE IF( BETA.NE.ONE )THEN 1901 DO 100, I = 1, J 1902 C( I, J ) = BETA*C( I, J ) 1903 100 CONTINUE 1904 END IF 1905 DO 120, L = 1, K 1906 IF( ( A( J, L ).NE.ZERO ).OR. 1907 $ ( B( J, L ).NE.ZERO ) )THEN 1908 TEMP1 = ALPHA*B( J, L ) 1909 TEMP2 = ALPHA*A( J, L ) 1910 DO 110, I = 1, J 1911 C( I, J ) = C( I, J ) + 1912 $ A( I, L )*TEMP1 + B( I, L )*TEMP2 1913 110 CONTINUE 1914 END IF 1915 120 CONTINUE 1916 130 CONTINUE 1917 ELSE 1918 DO 180, J = 1, N 1919 IF( BETA.EQ.ZERO )THEN 1920 DO 140, I = J, N 1921 C( I, J ) = ZERO 1922 140 CONTINUE 1923 ELSE IF( BETA.NE.ONE )THEN 1924 DO 150, I = J, N 1925 C( I, J ) = BETA*C( I, J ) 1926 150 CONTINUE 1927 END IF 1928 DO 170, L = 1, K 1929 IF( ( A( J, L ).NE.ZERO ).OR. 1930 $ ( B( J, L ).NE.ZERO ) )THEN 1931 TEMP1 = ALPHA*B( J, L ) 1932 TEMP2 = ALPHA*A( J, L ) 1933 DO 160, I = J, N 1934 C( I, J ) = C( I, J ) + 1935 $ A( I, L )*TEMP1 + B( I, L )*TEMP2 1936 160 CONTINUE 1937 END IF 1938 170 CONTINUE 1939 180 CONTINUE 1940 END IF 1941 ELSE 1942* 1943* Form C := alpha*A'*B + alpha*B'*A + C. 1944* 1945 IF( UPPER )THEN 1946 DO 210, J = 1, N 1947 DO 200, I = 1, J 1948 TEMP1 = ZERO 1949 TEMP2 = ZERO 1950 DO 190, L = 1, K 1951 TEMP1 = TEMP1 + A( L, I )*B( L, J ) 1952 TEMP2 = TEMP2 + B( L, I )*A( L, J ) 1953 190 CONTINUE 1954 IF( BETA.EQ.ZERO )THEN 1955 C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 1956 ELSE 1957 C( I, J ) = BETA *C( I, J ) + 1958 $ ALPHA*TEMP1 + ALPHA*TEMP2 1959 END IF 1960 200 CONTINUE 1961 210 CONTINUE 1962 ELSE 1963 DO 240, J = 1, N 1964 DO 230, I = J, N 1965 TEMP1 = ZERO 1966 TEMP2 = ZERO 1967 DO 220, L = 1, K 1968 TEMP1 = TEMP1 + A( L, I )*B( L, J ) 1969 TEMP2 = TEMP2 + B( L, I )*A( L, J ) 1970 220 CONTINUE 1971 IF( BETA.EQ.ZERO )THEN 1972 C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 1973 ELSE 1974 C( I, J ) = BETA *C( I, J ) + 1975 $ ALPHA*TEMP1 + ALPHA*TEMP2 1976 END IF 1977 230 CONTINUE 1978 240 CONTINUE 1979 END IF 1980 END IF 1981* 1982 RETURN 1983* 1984* End of DSYR2K. 1985* 1986 END 1987 SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, 1988 $ B, LDB ) 1989* .. Scalar Arguments .. 1990 CHARACTER*1 SIDE, UPLO, TRANSA, DIAG 1991 INTEGER M, N, LDA, LDB 1992 DOUBLE PRECISION ALPHA 1993* .. Array Arguments .. 1994 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 1995* .. 1996* 1997* Purpose 1998* ======= 1999* 2000* DTRMM performs one of the matrix-matrix operations 2001* 2002* B := alpha*op( A )*B, or B := alpha*B*op( A ), 2003* 2004* where alpha is a scalar, B is an m by n matrix, A is a unit, or 2005* non-unit, upper or lower triangular matrix and op( A ) is one of 2006* 2007* op( A ) = A or op( A ) = A'. 2008* 2009* Parameters 2010* ========== 2011* 2012* SIDE - CHARACTER*1. 2013* On entry, SIDE specifies whether op( A ) multiplies B from 2014* the left or right as follows: 2015* 2016* SIDE = 'L' or 'l' B := alpha*op( A )*B. 2017* 2018* SIDE = 'R' or 'r' B := alpha*B*op( A ). 2019* 2020* Unchanged on exit. 2021* 2022* UPLO - CHARACTER*1. 2023* On entry, UPLO specifies whether the matrix A is an upper or 2024* lower triangular matrix as follows: 2025* 2026* UPLO = 'U' or 'u' A is an upper triangular matrix. 2027* 2028* UPLO = 'L' or 'l' A is a lower triangular matrix. 2029* 2030* Unchanged on exit. 2031* 2032* TRANSA - CHARACTER*1. 2033* On entry, TRANSA specifies the form of op( A ) to be used in 2034* the matrix multiplication as follows: 2035* 2036* TRANSA = 'N' or 'n' op( A ) = A. 2037* 2038* TRANSA = 'T' or 't' op( A ) = A'. 2039* 2040* TRANSA = 'C' or 'c' op( A ) = A'. 2041* 2042* Unchanged on exit. 2043* 2044* DIAG - CHARACTER*1. 2045* On entry, DIAG specifies whether or not A is unit triangular 2046* as follows: 2047* 2048* DIAG = 'U' or 'u' A is assumed to be unit triangular. 2049* 2050* DIAG = 'N' or 'n' A is not assumed to be unit 2051* triangular. 2052* 2053* Unchanged on exit. 2054* 2055* M - INTEGER. 2056* On entry, M specifies the number of rows of B. M must be at 2057* least zero. 2058* Unchanged on exit. 2059* 2060* N - INTEGER. 2061* On entry, N specifies the number of columns of B. N must be 2062* at least zero. 2063* Unchanged on exit. 2064* 2065* ALPHA - DOUBLE PRECISION. 2066* On entry, ALPHA specifies the scalar alpha. When alpha is 2067* zero then A is not referenced and B need not be set before 2068* entry. 2069* Unchanged on exit. 2070* 2071* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m 2072* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. 2073* Before entry with UPLO = 'U' or 'u', the leading k by k 2074* upper triangular part of the array A must contain the upper 2075* triangular matrix and the strictly lower triangular part of 2076* A is not referenced. 2077* Before entry with UPLO = 'L' or 'l', the leading k by k 2078* lower triangular part of the array A must contain the lower 2079* triangular matrix and the strictly upper triangular part of 2080* A is not referenced. 2081* Note that when DIAG = 'U' or 'u', the diagonal elements of 2082* A are not referenced either, but are assumed to be unity. 2083* Unchanged on exit. 2084* 2085* LDA - INTEGER. 2086* On entry, LDA specifies the first dimension of A as declared 2087* in the calling (sub) program. When SIDE = 'L' or 'l' then 2088* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' 2089* then LDA must be at least max( 1, n ). 2090* Unchanged on exit. 2091* 2092* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). 2093* Before entry, the leading m by n part of the array B must 2094* contain the matrix B, and on exit is overwritten by the 2095* transformed matrix. 2096* 2097* LDB - INTEGER. 2098* On entry, LDB specifies the first dimension of B as declared 2099* in the calling (sub) program. LDB must be at least 2100* max( 1, m ). 2101* Unchanged on exit. 2102* 2103* 2104* Level 3 Blas routine. 2105* 2106* -- Written on 8-February-1989. 2107* Jack Dongarra, Argonne National Laboratory. 2108* Iain Duff, AERE Harwell. 2109* Jeremy Du Croz, Numerical Algorithms Group Ltd. 2110* Sven Hammarling, Numerical Algorithms Group Ltd. 2111* 2112* 2113* .. External Functions .. 2114 LOGICAL LSAME 2115 EXTERNAL LSAME 2116* .. External Subroutines .. 2117 EXTERNAL XERBLA 2118* .. Intrinsic Functions .. 2119 INTRINSIC MAX 2120* .. Local Scalars .. 2121 LOGICAL LSIDE, NOUNIT, UPPER 2122 INTEGER I, INFO, J, K, NROWA 2123 DOUBLE PRECISION TEMP 2124* .. Parameters .. 2125 DOUBLE PRECISION ONE , ZERO 2126 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 2127* .. 2128* .. Executable Statements .. 2129* 2130* Test the input parameters. 2131* 2132 LSIDE = LSAME( SIDE , 'L' ) 2133 IF( LSIDE )THEN 2134 NROWA = M 2135 ELSE 2136 NROWA = N 2137 END IF 2138 NOUNIT = LSAME( DIAG , 'N' ) 2139 UPPER = LSAME( UPLO , 'U' ) 2140* 2141 INFO = 0 2142 IF( ( .NOT.LSIDE ).AND. 2143 $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN 2144 INFO = 1 2145 ELSE IF( ( .NOT.UPPER ).AND. 2146 $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN 2147 INFO = 2 2148 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. 2149 $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. 2150 $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN 2151 INFO = 3 2152 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. 2153 $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN 2154 INFO = 4 2155 ELSE IF( M .LT.0 )THEN 2156 INFO = 5 2157 ELSE IF( N .LT.0 )THEN 2158 INFO = 6 2159 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN 2160 INFO = 9 2161 ELSE IF( LDB.LT.MAX( 1, M ) )THEN 2162 INFO = 11 2163 END IF 2164 IF( INFO.NE.0 )THEN 2165 CALL XERBLA( 'DTRMM ', INFO ) 2166 RETURN 2167 END IF 2168* 2169* Quick return if possible. 2170* 2171 IF( N.EQ.0 ) 2172 $ RETURN 2173* 2174* And when alpha.eq.zero. 2175* 2176 IF( ALPHA.EQ.ZERO )THEN 2177 DO 20, J = 1, N 2178 DO 10, I = 1, M 2179 B( I, J ) = ZERO 2180 10 CONTINUE 2181 20 CONTINUE 2182 RETURN 2183 END IF 2184* 2185* Start the operations. 2186* 2187 IF( LSIDE )THEN 2188 IF( LSAME( TRANSA, 'N' ) )THEN 2189* 2190* Form B := alpha*A*B. 2191* 2192 IF( UPPER )THEN 2193 DO 50, J = 1, N 2194 DO 40, K = 1, M 2195 IF( B( K, J ).NE.ZERO )THEN 2196 TEMP = ALPHA*B( K, J ) 2197 DO 30, I = 1, K - 1 2198 B( I, J ) = B( I, J ) + TEMP*A( I, K ) 2199 30 CONTINUE 2200 IF( NOUNIT ) 2201 $ TEMP = TEMP*A( K, K ) 2202 B( K, J ) = TEMP 2203 END IF 2204 40 CONTINUE 2205 50 CONTINUE 2206 ELSE 2207 DO 80, J = 1, N 2208 DO 70 K = M, 1, -1 2209 IF( B( K, J ).NE.ZERO )THEN 2210 TEMP = ALPHA*B( K, J ) 2211 B( K, J ) = TEMP 2212 IF( NOUNIT ) 2213 $ B( K, J ) = B( K, J )*A( K, K ) 2214 DO 60, I = K + 1, M 2215 B( I, J ) = B( I, J ) + TEMP*A( I, K ) 2216 60 CONTINUE 2217 END IF 2218 70 CONTINUE 2219 80 CONTINUE 2220 END IF 2221 ELSE 2222* 2223* Form B := alpha*B*A'. 2224* 2225 IF( UPPER )THEN 2226 DO 110, J = 1, N 2227 DO 100, I = M, 1, -1 2228 TEMP = B( I, J ) 2229 IF( NOUNIT ) 2230 $ TEMP = TEMP*A( I, I ) 2231 DO 90, K = 1, I - 1 2232 TEMP = TEMP + A( K, I )*B( K, J ) 2233 90 CONTINUE 2234 B( I, J ) = ALPHA*TEMP 2235 100 CONTINUE 2236 110 CONTINUE 2237 ELSE 2238 DO 140, J = 1, N 2239 DO 130, I = 1, M 2240 TEMP = B( I, J ) 2241 IF( NOUNIT ) 2242 $ TEMP = TEMP*A( I, I ) 2243 DO 120, K = I + 1, M 2244 TEMP = TEMP + A( K, I )*B( K, J ) 2245 120 CONTINUE 2246 B( I, J ) = ALPHA*TEMP 2247 130 CONTINUE 2248 140 CONTINUE 2249 END IF 2250 END IF 2251 ELSE 2252 IF( LSAME( TRANSA, 'N' ) )THEN 2253* 2254* Form B := alpha*B*A. 2255* 2256 IF( UPPER )THEN 2257 DO 180, J = N, 1, -1 2258 TEMP = ALPHA 2259 IF( NOUNIT ) 2260 $ TEMP = TEMP*A( J, J ) 2261 DO 150, I = 1, M 2262 B( I, J ) = TEMP*B( I, J ) 2263 150 CONTINUE 2264 DO 170, K = 1, J - 1 2265 IF( A( K, J ).NE.ZERO )THEN 2266 TEMP = ALPHA*A( K, J ) 2267 DO 160, I = 1, M 2268 B( I, J ) = B( I, J ) + TEMP*B( I, K ) 2269 160 CONTINUE 2270 END IF 2271 170 CONTINUE 2272 180 CONTINUE 2273 ELSE 2274 DO 220, J = 1, N 2275 TEMP = ALPHA 2276 IF( NOUNIT ) 2277 $ TEMP = TEMP*A( J, J ) 2278 DO 190, I = 1, M 2279 B( I, J ) = TEMP*B( I, J ) 2280 190 CONTINUE 2281 DO 210, K = J + 1, N 2282 IF( A( K, J ).NE.ZERO )THEN 2283 TEMP = ALPHA*A( K, J ) 2284 DO 200, I = 1, M 2285 B( I, J ) = B( I, J ) + TEMP*B( I, K ) 2286 200 CONTINUE 2287 END IF 2288 210 CONTINUE 2289 220 CONTINUE 2290 END IF 2291 ELSE 2292* 2293* Form B := alpha*B*A'. 2294* 2295 IF( UPPER )THEN 2296 DO 260, K = 1, N 2297 DO 240, J = 1, K - 1 2298 IF( A( J, K ).NE.ZERO )THEN 2299 TEMP = ALPHA*A( J, K ) 2300 DO 230, I = 1, M 2301 B( I, J ) = B( I, J ) + TEMP*B( I, K ) 2302 230 CONTINUE 2303 END IF 2304 240 CONTINUE 2305 TEMP = ALPHA 2306 IF( NOUNIT ) 2307 $ TEMP = TEMP*A( K, K ) 2308 IF( TEMP.NE.ONE )THEN 2309 DO 250, I = 1, M 2310 B( I, K ) = TEMP*B( I, K ) 2311 250 CONTINUE 2312 END IF 2313 260 CONTINUE 2314 ELSE 2315 DO 300, K = N, 1, -1 2316 DO 280, J = K + 1, N 2317 IF( A( J, K ).NE.ZERO )THEN 2318 TEMP = ALPHA*A( J, K ) 2319 DO 270, I = 1, M 2320 B( I, J ) = B( I, J ) + TEMP*B( I, K ) 2321 270 CONTINUE 2322 END IF 2323 280 CONTINUE 2324 TEMP = ALPHA 2325 IF( NOUNIT ) 2326 $ TEMP = TEMP*A( K, K ) 2327 IF( TEMP.NE.ONE )THEN 2328 DO 290, I = 1, M 2329 B( I, K ) = TEMP*B( I, K ) 2330 290 CONTINUE 2331 END IF 2332 300 CONTINUE 2333 END IF 2334 END IF 2335 END IF 2336* 2337 RETURN 2338* 2339* End of DTRMM . 2340* 2341 END 2342 SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) 2343* .. Scalar Arguments .. 2344 INTEGER INCX, LDA, N 2345 CHARACTER*1 DIAG, TRANS, UPLO 2346* .. Array Arguments .. 2347 DOUBLE PRECISION A( LDA, * ), X( * ) 2348* .. 2349* 2350* Purpose 2351* ======= 2352* 2353* DTRMV performs one of the matrix-vector operations 2354* 2355* x := A*x, or x := A'*x, 2356* 2357* where x is an n element vector and A is an n by n unit, or non-unit, 2358* upper or lower triangular matrix. 2359* 2360* Parameters 2361* ========== 2362* 2363* UPLO - CHARACTER*1. 2364* On entry, UPLO specifies whether the matrix is an upper or 2365* lower triangular matrix as follows: 2366* 2367* UPLO = 'U' or 'u' A is an upper triangular matrix. 2368* 2369* UPLO = 'L' or 'l' A is a lower triangular matrix. 2370* 2371* Unchanged on exit. 2372* 2373* TRANS - CHARACTER*1. 2374* On entry, TRANS specifies the operation to be performed as 2375* follows: 2376* 2377* TRANS = 'N' or 'n' x := A*x. 2378* 2379* TRANS = 'T' or 't' x := A'*x. 2380* 2381* TRANS = 'C' or 'c' x := A'*x. 2382* 2383* Unchanged on exit. 2384* 2385* DIAG - CHARACTER*1. 2386* On entry, DIAG specifies whether or not A is unit 2387* triangular as follows: 2388* 2389* DIAG = 'U' or 'u' A is assumed to be unit triangular. 2390* 2391* DIAG = 'N' or 'n' A is not assumed to be unit 2392* triangular. 2393* 2394* Unchanged on exit. 2395* 2396* N - INTEGER. 2397* On entry, N specifies the order of the matrix A. 2398* N must be at least zero. 2399* Unchanged on exit. 2400* 2401* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 2402* Before entry with UPLO = 'U' or 'u', the leading n by n 2403* upper triangular part of the array A must contain the upper 2404* triangular matrix and the strictly lower triangular part of 2405* A is not referenced. 2406* Before entry with UPLO = 'L' or 'l', the leading n by n 2407* lower triangular part of the array A must contain the lower 2408* triangular matrix and the strictly upper triangular part of 2409* A is not referenced. 2410* Note that when DIAG = 'U' or 'u', the diagonal elements of 2411* A are not referenced either, but are assumed to be unity. 2412* Unchanged on exit. 2413* 2414* LDA - INTEGER. 2415* On entry, LDA specifies the first dimension of A as declared 2416* in the calling (sub) program. LDA must be at least 2417* max( 1, n ). 2418* Unchanged on exit. 2419* 2420* X - DOUBLE PRECISION array of dimension at least 2421* ( 1 + ( n - 1 )*abs( INCX ) ). 2422* Before entry, the incremented array X must contain the n 2423* element vector x. On exit, X is overwritten with the 2424* tranformed vector x. 2425* 2426* INCX - INTEGER. 2427* On entry, INCX specifies the increment for the elements of 2428* X. INCX must not be zero. 2429* Unchanged on exit. 2430* 2431* 2432* Level 2 Blas routine. 2433* 2434* -- Written on 22-October-1986. 2435* Jack Dongarra, Argonne National Lab. 2436* Jeremy Du Croz, Nag Central Office. 2437* Sven Hammarling, Nag Central Office. 2438* Richard Hanson, Sandia National Labs. 2439* 2440* 2441* .. Parameters .. 2442 DOUBLE PRECISION ZERO 2443 PARAMETER ( ZERO = 0.0D+0 ) 2444* .. Local Scalars .. 2445 DOUBLE PRECISION TEMP 2446 INTEGER I, INFO, IX, J, JX, KX 2447 LOGICAL NOUNIT 2448* .. External Functions .. 2449 LOGICAL LSAME 2450 EXTERNAL LSAME 2451* .. External Subroutines .. 2452 EXTERNAL XERBLA 2453* .. Intrinsic Functions .. 2454 INTRINSIC MAX 2455* .. 2456* .. Executable Statements .. 2457* 2458* Test the input parameters. 2459* 2460 INFO = 0 2461 IF ( .NOT.LSAME( UPLO , 'U' ).AND. 2462 $ .NOT.LSAME( UPLO , 'L' ) )THEN 2463 INFO = 1 2464 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. 2465 $ .NOT.LSAME( TRANS, 'T' ).AND. 2466 $ .NOT.LSAME( TRANS, 'C' ) )THEN 2467 INFO = 2 2468 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. 2469 $ .NOT.LSAME( DIAG , 'N' ) )THEN 2470 INFO = 3 2471 ELSE IF( N.LT.0 )THEN 2472 INFO = 4 2473 ELSE IF( LDA.LT.MAX( 1, N ) )THEN 2474 INFO = 6 2475 ELSE IF( INCX.EQ.0 )THEN 2476 INFO = 8 2477 END IF 2478 IF( INFO.NE.0 )THEN 2479 CALL XERBLA( 'DTRMV ', INFO ) 2480 RETURN 2481 END IF 2482* 2483* Quick return if possible. 2484* 2485 IF( N.EQ.0 ) 2486 $ RETURN 2487* 2488 NOUNIT = LSAME( DIAG, 'N' ) 2489* 2490* Set up the start point in X if the increment is not unity. This 2491* will be ( N - 1 )*INCX too small for descending loops. 2492* 2493 IF( INCX.LE.0 )THEN 2494 KX = 1 - ( N - 1 )*INCX 2495 ELSE IF( INCX.NE.1 )THEN 2496 KX = 1 2497 END IF 2498* 2499* Start the operations. In this version the elements of A are 2500* accessed sequentially with one pass through A. 2501* 2502 IF( LSAME( TRANS, 'N' ) )THEN 2503* 2504* Form x := A*x. 2505* 2506 IF( LSAME( UPLO, 'U' ) )THEN 2507 IF( INCX.EQ.1 )THEN 2508 DO 20, J = 1, N 2509 IF( X( J ).NE.ZERO )THEN 2510 TEMP = X( J ) 2511 DO 10, I = 1, J - 1 2512 X( I ) = X( I ) + TEMP*A( I, J ) 2513 10 CONTINUE 2514 IF( NOUNIT ) 2515 $ X( J ) = X( J )*A( J, J ) 2516 END IF 2517 20 CONTINUE 2518 ELSE 2519 JX = KX 2520 DO 40, J = 1, N 2521 IF( X( JX ).NE.ZERO )THEN 2522 TEMP = X( JX ) 2523 IX = KX 2524 DO 30, I = 1, J - 1 2525 X( IX ) = X( IX ) + TEMP*A( I, J ) 2526 IX = IX + INCX 2527 30 CONTINUE 2528 IF( NOUNIT ) 2529 $ X( JX ) = X( JX )*A( J, J ) 2530 END IF 2531 JX = JX + INCX 2532 40 CONTINUE 2533 END IF 2534 ELSE 2535 IF( INCX.EQ.1 )THEN 2536 DO 60, J = N, 1, -1 2537 IF( X( J ).NE.ZERO )THEN 2538 TEMP = X( J ) 2539 DO 50, I = N, J + 1, -1 2540 X( I ) = X( I ) + TEMP*A( I, J ) 2541 50 CONTINUE 2542 IF( NOUNIT ) 2543 $ X( J ) = X( J )*A( J, J ) 2544 END IF 2545 60 CONTINUE 2546 ELSE 2547 KX = KX + ( N - 1 )*INCX 2548 JX = KX 2549 DO 80, J = N, 1, -1 2550 IF( X( JX ).NE.ZERO )THEN 2551 TEMP = X( JX ) 2552 IX = KX 2553 DO 70, I = N, J + 1, -1 2554 X( IX ) = X( IX ) + TEMP*A( I, J ) 2555 IX = IX - INCX 2556 70 CONTINUE 2557 IF( NOUNIT ) 2558 $ X( JX ) = X( JX )*A( J, J ) 2559 END IF 2560 JX = JX - INCX 2561 80 CONTINUE 2562 END IF 2563 END IF 2564 ELSE 2565* 2566* Form x := A'*x. 2567* 2568 IF( LSAME( UPLO, 'U' ) )THEN 2569 IF( INCX.EQ.1 )THEN 2570 DO 100, J = N, 1, -1 2571 TEMP = X( J ) 2572 IF( NOUNIT ) 2573 $ TEMP = TEMP*A( J, J ) 2574 DO 90, I = J - 1, 1, -1 2575 TEMP = TEMP + A( I, J )*X( I ) 2576 90 CONTINUE 2577 X( J ) = TEMP 2578 100 CONTINUE 2579 ELSE 2580 JX = KX + ( N - 1 )*INCX 2581 DO 120, J = N, 1, -1 2582 TEMP = X( JX ) 2583 IX = JX 2584 IF( NOUNIT ) 2585 $ TEMP = TEMP*A( J, J ) 2586 DO 110, I = J - 1, 1, -1 2587 IX = IX - INCX 2588 TEMP = TEMP + A( I, J )*X( IX ) 2589 110 CONTINUE 2590 X( JX ) = TEMP 2591 JX = JX - INCX 2592 120 CONTINUE 2593 END IF 2594 ELSE 2595 IF( INCX.EQ.1 )THEN 2596 DO 140, J = 1, N 2597 TEMP = X( J ) 2598 IF( NOUNIT ) 2599 $ TEMP = TEMP*A( J, J ) 2600 DO 130, I = J + 1, N 2601 TEMP = TEMP + A( I, J )*X( I ) 2602 130 CONTINUE 2603 X( J ) = TEMP 2604 140 CONTINUE 2605 ELSE 2606 JX = KX 2607 DO 160, J = 1, N 2608 TEMP = X( JX ) 2609 IX = JX 2610 IF( NOUNIT ) 2611 $ TEMP = TEMP*A( J, J ) 2612 DO 150, I = J + 1, N 2613 IX = IX + INCX 2614 TEMP = TEMP + A( I, J )*X( IX ) 2615 150 CONTINUE 2616 X( JX ) = TEMP 2617 JX = JX + INCX 2618 160 CONTINUE 2619 END IF 2620 END IF 2621 END IF 2622* 2623 RETURN 2624* 2625* End of DTRMV . 2626* 2627 END 2628 SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) 2629* .. Scalar Arguments .. 2630 INTEGER INCX, LDA, N 2631 CHARACTER*1 DIAG, TRANS, UPLO 2632* .. Array Arguments .. 2633 DOUBLE PRECISION A( LDA, * ), X( * ) 2634* .. 2635* 2636* Purpose 2637* ======= 2638* 2639* DTRSV solves one of the systems of equations 2640* 2641* A*x = b, or A'*x = b, 2642* 2643* where b and x are n element vectors and A is an n by n unit, or 2644* non-unit, upper or lower triangular matrix. 2645* 2646* No test for singularity or near-singularity is included in this 2647* routine. Such tests must be performed before calling this routine. 2648* 2649* Parameters 2650* ========== 2651* 2652* UPLO - CHARACTER*1. 2653* On entry, UPLO specifies whether the matrix is an upper or 2654* lower triangular matrix as follows: 2655* 2656* UPLO = 'U' or 'u' A is an upper triangular matrix. 2657* 2658* UPLO = 'L' or 'l' A is a lower triangular matrix. 2659* 2660* Unchanged on exit. 2661* 2662* TRANS - CHARACTER*1. 2663* On entry, TRANS specifies the equations to be solved as 2664* follows: 2665* 2666* TRANS = 'N' or 'n' A*x = b. 2667* 2668* TRANS = 'T' or 't' A'*x = b. 2669* 2670* TRANS = 'C' or 'c' A'*x = b. 2671* 2672* Unchanged on exit. 2673* 2674* DIAG - CHARACTER*1. 2675* On entry, DIAG specifies whether or not A is unit 2676* triangular as follows: 2677* 2678* DIAG = 'U' or 'u' A is assumed to be unit triangular. 2679* 2680* DIAG = 'N' or 'n' A is not assumed to be unit 2681* triangular. 2682* 2683* Unchanged on exit. 2684* 2685* N - INTEGER. 2686* On entry, N specifies the order of the matrix A. 2687* N must be at least zero. 2688* Unchanged on exit. 2689* 2690* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 2691* Before entry with UPLO = 'U' or 'u', the leading n by n 2692* upper triangular part of the array A must contain the upper 2693* triangular matrix and the strictly lower triangular part of 2694* A is not referenced. 2695* Before entry with UPLO = 'L' or 'l', the leading n by n 2696* lower triangular part of the array A must contain the lower 2697* triangular matrix and the strictly upper triangular part of 2698* A is not referenced. 2699* Note that when DIAG = 'U' or 'u', the diagonal elements of 2700* A are not referenced either, but are assumed to be unity. 2701* Unchanged on exit. 2702* 2703* LDA - INTEGER. 2704* On entry, LDA specifies the first dimension of A as declared 2705* in the calling (sub) program. LDA must be at least 2706* max( 1, n ). 2707* Unchanged on exit. 2708* 2709* X - DOUBLE PRECISION array of dimension at least 2710* ( 1 + ( n - 1 )*abs( INCX ) ). 2711* Before entry, the incremented array X must contain the n 2712* element right-hand side vector b. On exit, X is overwritten 2713* with the solution vector x. 2714* 2715* INCX - INTEGER. 2716* On entry, INCX specifies the increment for the elements of 2717* X. INCX must not be zero. 2718* Unchanged on exit. 2719* 2720* 2721* Level 2 Blas routine. 2722* 2723* -- Written on 22-October-1986. 2724* Jack Dongarra, Argonne National Lab. 2725* Jeremy Du Croz, Nag Central Office. 2726* Sven Hammarling, Nag Central Office. 2727* Richard Hanson, Sandia National Labs. 2728* 2729* 2730* .. Parameters .. 2731 DOUBLE PRECISION ZERO 2732 PARAMETER ( ZERO = 0.0D+0 ) 2733* .. Local Scalars .. 2734 DOUBLE PRECISION TEMP 2735 INTEGER I, INFO, IX, J, JX, KX 2736 LOGICAL NOUNIT 2737* .. External Functions .. 2738 LOGICAL LSAME 2739 EXTERNAL LSAME 2740* .. External Subroutines .. 2741 EXTERNAL XERBLA 2742* .. Intrinsic Functions .. 2743 INTRINSIC MAX 2744* .. 2745* .. Executable Statements .. 2746* 2747* Test the input parameters. 2748* 2749 INFO = 0 2750 IF ( .NOT.LSAME( UPLO , 'U' ).AND. 2751 $ .NOT.LSAME( UPLO , 'L' ) )THEN 2752 INFO = 1 2753 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. 2754 $ .NOT.LSAME( TRANS, 'T' ).AND. 2755 $ .NOT.LSAME( TRANS, 'C' ) )THEN 2756 INFO = 2 2757 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. 2758 $ .NOT.LSAME( DIAG , 'N' ) )THEN 2759 INFO = 3 2760 ELSE IF( N.LT.0 )THEN 2761 INFO = 4 2762 ELSE IF( LDA.LT.MAX( 1, N ) )THEN 2763 INFO = 6 2764 ELSE IF( INCX.EQ.0 )THEN 2765 INFO = 8 2766 END IF 2767 IF( INFO.NE.0 )THEN 2768 CALL XERBLA( 'DTRSV ', INFO ) 2769 RETURN 2770 END IF 2771* 2772* Quick return if possible. 2773* 2774 IF( N.EQ.0 ) 2775 $ RETURN 2776* 2777 NOUNIT = LSAME( DIAG, 'N' ) 2778* 2779* Set up the start point in X if the increment is not unity. This 2780* will be ( N - 1 )*INCX too small for descending loops. 2781* 2782 IF( INCX.LE.0 )THEN 2783 KX = 1 - ( N - 1 )*INCX 2784 ELSE IF( INCX.NE.1 )THEN 2785 KX = 1 2786 END IF 2787* 2788* Start the operations. In this version the elements of A are 2789* accessed sequentially with one pass through A. 2790* 2791 IF( LSAME( TRANS, 'N' ) )THEN 2792* 2793* Form x := inv( A )*x. 2794* 2795 IF( LSAME( UPLO, 'U' ) )THEN 2796 IF( INCX.EQ.1 )THEN 2797 DO 20, J = N, 1, -1 2798 IF( X( J ).NE.ZERO )THEN 2799 IF( NOUNIT ) 2800 $ X( J ) = X( J )/A( J, J ) 2801 TEMP = X( J ) 2802 DO 10, I = J - 1, 1, -1 2803 X( I ) = X( I ) - TEMP*A( I, J ) 2804 10 CONTINUE 2805 END IF 2806 20 CONTINUE 2807 ELSE 2808 JX = KX + ( N - 1 )*INCX 2809 DO 40, J = N, 1, -1 2810 IF( X( JX ).NE.ZERO )THEN 2811 IF( NOUNIT ) 2812 $ X( JX ) = X( JX )/A( J, J ) 2813 TEMP = X( JX ) 2814 IX = JX 2815 DO 30, I = J - 1, 1, -1 2816 IX = IX - INCX 2817 X( IX ) = X( IX ) - TEMP*A( I, J ) 2818 30 CONTINUE 2819 END IF 2820 JX = JX - INCX 2821 40 CONTINUE 2822 END IF 2823 ELSE 2824 IF( INCX.EQ.1 )THEN 2825 DO 60, J = 1, N 2826 IF( X( J ).NE.ZERO )THEN 2827 IF( NOUNIT ) 2828 $ X( J ) = X( J )/A( J, J ) 2829 TEMP = X( J ) 2830 DO 50, I = J + 1, N 2831 X( I ) = X( I ) - TEMP*A( I, J ) 2832 50 CONTINUE 2833 END IF 2834 60 CONTINUE 2835 ELSE 2836 JX = KX 2837 DO 80, J = 1, N 2838 IF( X( JX ).NE.ZERO )THEN 2839 IF( NOUNIT ) 2840 $ X( JX ) = X( JX )/A( J, J ) 2841 TEMP = X( JX ) 2842 IX = JX 2843 DO 70, I = J + 1, N 2844 IX = IX + INCX 2845 X( IX ) = X( IX ) - TEMP*A( I, J ) 2846 70 CONTINUE 2847 END IF 2848 JX = JX + INCX 2849 80 CONTINUE 2850 END IF 2851 END IF 2852 ELSE 2853* 2854* Form x := inv( A' )*x. 2855* 2856 IF( LSAME( UPLO, 'U' ) )THEN 2857 IF( INCX.EQ.1 )THEN 2858 DO 100, J = 1, N 2859 TEMP = X( J ) 2860 DO 90, I = 1, J - 1 2861 TEMP = TEMP - A( I, J )*X( I ) 2862 90 CONTINUE 2863 IF( NOUNIT ) 2864 $ TEMP = TEMP/A( J, J ) 2865 X( J ) = TEMP 2866 100 CONTINUE 2867 ELSE 2868 JX = KX 2869 DO 120, J = 1, N 2870 TEMP = X( JX ) 2871 IX = KX 2872 DO 110, I = 1, J - 1 2873 TEMP = TEMP - A( I, J )*X( IX ) 2874 IX = IX + INCX 2875 110 CONTINUE 2876 IF( NOUNIT ) 2877 $ TEMP = TEMP/A( J, J ) 2878 X( JX ) = TEMP 2879 JX = JX + INCX 2880 120 CONTINUE 2881 END IF 2882 ELSE 2883 IF( INCX.EQ.1 )THEN 2884 DO 140, J = N, 1, -1 2885 TEMP = X( J ) 2886 DO 130, I = N, J + 1, -1 2887 TEMP = TEMP - A( I, J )*X( I ) 2888 130 CONTINUE 2889 IF( NOUNIT ) 2890 $ TEMP = TEMP/A( J, J ) 2891 X( J ) = TEMP 2892 140 CONTINUE 2893 ELSE 2894 KX = KX + ( N - 1 )*INCX 2895 JX = KX 2896 DO 160, J = N, 1, -1 2897 TEMP = X( JX ) 2898 IX = KX 2899 DO 150, I = N, J + 1, -1 2900 TEMP = TEMP - A( I, J )*X( IX ) 2901 IX = IX - INCX 2902 150 CONTINUE 2903 IF( NOUNIT ) 2904 $ TEMP = TEMP/A( J, J ) 2905 X( JX ) = TEMP 2906 JX = JX - INCX 2907 160 CONTINUE 2908 END IF 2909 END IF 2910 END IF 2911* 2912 RETURN 2913* 2914* End of DTRSV . 2915* 2916 END 2917 integer function idamax(n,dx,incx) 2918c 2919c finds the index of element having max. absolute value. 2920c jack dongarra, linpack, 3/11/78. 2921c modified 3/93 to return if incx .le. 0. 2922c 2923 double precision dx(1),dmax 2924 integer i,incx,ix,n 2925c 2926 idamax = 0 2927 if( n.lt.1 .or. incx.le.0 ) return 2928 idamax = 1 2929 if(n.eq.1)return 2930 if(incx.eq.1)go to 20 2931c 2932c code for increment not equal to 1 2933c 2934 ix = 1 2935 dmax = dabs(dx(1)) 2936 ix = ix + incx 2937 do 10 i = 2,n 2938 if(dabs(dx(ix)).le.dmax) go to 5 2939 idamax = i 2940 dmax = dabs(dx(ix)) 2941 5 ix = ix + incx 2942 10 continue 2943 return 2944c 2945c code for increment equal to 1 2946c 2947 20 dmax = dabs(dx(1)) 2948 do 30 i = 2,n 2949 if(dabs(dx(i)).le.dmax) go to 30 2950 idamax = i 2951 dmax = dabs(dx(i)) 2952 30 continue 2953 return 2954 end 2955