1C 2C 3C... Dalton, Release DALTON2016 4C... 5C... These routines are in the public domain and can be 6C... used freely in other programs. 7C... 8C 9C FILE : pdpack/linextra.F 10C 11C Extra matrix and vector manipulation routines and driver routines, 12C which not (always) are included in blas and linpack. 13C 14C Most of them are written by Hans Jørgen Aa. Jensen in the late 1980's. 15C These routines may be freely used by anyone. 16C 17C 1) DNORM2 (emulate ESSL DNORM2: do not use extended precision for intermediates) 18C 2) DAPTGE, DGETAP etc. for triangular packing and unpacking of matrices 19C (AP: Antisymmetric Packed, SP, Symmetric Packed, etc. 20C (ex: DAPTGE: Antisymmetric Packed To GEneral matrix) 21C NOTE: packed part is always LOWER triangle (important for AP) 22C 3) DSUM, IDAMIN, IDMAX, IDMIN (supplement DASUM and IDAMAX) 23C 4) "iblas": ICOPY, ISWAP, ... 24C 5) other extra routines: DGEZERO, DSPZERO, NDXGTA 25C (zero block of matrix; find next element .gt. a) 26C 6) DSPSOL, DSPSLI, DGESOL, DGEINV : driver routines calling LINPACK 27C 7) BLAS and LAPACK routines missing for ESSL: 28C LSAME; ILAENV, DSPTRI, DSPCON 29C 30C 31C 32C Block 1: 33C DNORM2 (emulate ESSL DNORM2: do not use extended precision for intermediates) 34C 35#if !defined (VAR_ESSL) && !defined (VAR_DXML) 36C DNORM2 is in ESSL library from IBM and DXML library on DEC-ALPHA 37C /* Deck dnorm2 */ 38 FUNCTION DNORM2(N,DX,INCX) 39C 40C Forms the two-norm of a vector. 41C 19-Sep-1988 -- hjaaj -- based on DNRM2 from LINPACK 42C This version does not use extended precision for intermediates 43C as the LINPACK version does.C 1) DNORM2 (emulate ESSL DNORM2: do not use extended precision for intermediates) 44 45C Equivalent to DNORM2 in IBM's ESSL library. 46C 47C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. 48C DNRM2: JACK DONGARRA, LINPACK, 3/11/78. 49C 50#include "implicit.h" 51C 52 DIMENSION DX(*) 53 PARAMETER ( ZERO = 0.0D0 ) 54C 55 IF (N.LE.0) THEN 56 DNORM2 = ZERO 57 RETURN 58 END IF 59 DTEMP = ZERO 60 IF(INCX.EQ.1)GO TO 20 61C 62C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 63C NOT EQUAL TO 1 64C 65 IX = 1 66 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 67 DO 10 I = 1,N 68 DTEMP = DTEMP + DX(IX)*DX(IX) 69 IX = IX + INCX 70 10 CONTINUE 71 DNORM2 = SQRT(DTEMP) 72 RETURN 73C 74C CODE FOR BOTH INCREMENTS EQUAL TO 1 75C 76C 77C CLEAN-UP LOOP 78C 79 20 M = MOD(N,5) 80 IF( M .EQ. 0 ) GO TO 40 81 DO 30 I = 1,M 82 DTEMP = DTEMP + DX(I)*DX(I) 83 30 CONTINUE 84 IF( N .LT. 5 ) GO TO 60 85 40 MP1 = M + 1 86 DO 50 I = MP1,N,5 87 DTEMP = DTEMP + DX(I)*DX(I) + DX(I + 1)*DX(I + 1) + 88 * DX(I + 2)*DX(I + 2) + DX(I + 3)*DX(I + 3) + DX(I + 4)*DX(I + 4) 89 50 CONTINUE 90 60 DNORM2 = SQRT(DTEMP) 91 RETURN 92 END 93#endif /* not VAR_ESSL */ 94C Block 2: 95C DAPTGE, DGETAP etc. for triangular packing and unpacking of matrices 96C (AP: Antisymmetric Packed, SP, Symmetric Packed, etc. 97C (ex: DAPTGE: Antisymmetric Packed To GEneral matrix) 98C NOTE: packed part is always LOWER triangle (important for AP) 99C /* Deck daptge */ 100 SUBROUTINE DAPTGE(N,AAP,AGE) 101C 102C 8-Feb-1987 Hans Joergen Aa. Jensen 103C 104C Purpose: Transform from AP format to GE format, that is: 105C unpack antisymmetric, packed (lower triangle) matrix AAP 106C to antisymmetric, unpacked matrix AGE. 107C 108#include "implicit.h" 109 DIMENSION AAP(*), AGE(N,*) 110C 111 DO 200 J = 1,N 112 JOFF = (J*J-J)/2 113 DO 100 I = 1,J-1 114 AGE(I,J) = - AAP(JOFF+I) 115 AGE(J,I) = AAP(JOFF+I) 116 100 CONTINUE 117 AGE(J,J) = AAP(JOFF+J) 118C ... is zero but included such that error may be detected. 119 200 CONTINUE 120C 121 RETURN 122 END 123C /* Deck dsptsi */ 124 SUBROUTINE DSPTSI(N,ASP,ASI) 125C 126C 8-Feb-1987 Hans Joergen Aa. Jensen 127C 128C Purpose: Transform from SP format to SI format, that is: 129C unpack symmetric, packed matrix ASP 130C to symmetric, unpacked matrix ASI. 131C 132#include "implicit.h" 133 DIMENSION ASP(*), ASI(N,*) 134C 135 ENTRY DSPTGE(N,ASP,ASI) 136C ... equivalent subroutine name 137C 138 DO 200 J = 1,N 139 JOFF = (J*J-J)/2 140 DO 100 I = 1,J-1 141 ASI(I,J) = ASP(JOFF+I) 142 ASI(J,I) = ASP(JOFF+I) 143 100 CONTINUE 144 ASI(J,J) = ASP(JOFF+J) 145 200 CONTINUE 146C 147 RETURN 148 END 149C /* Deck dgetap */ 150 SUBROUTINE DGETAP(N,AGE,AAP) 151C 152C 8-Feb-1987 Hans Joergen Aa. Jensen 153Cdgetap 154C Purpose: Transform from GE format to AP format, that is: 155C extract antisymmetric part of general matrix AGE 156C to antisymmetric, packed matrix AAP (lower 157C triangle saved). 158C 159#include "implicit.h" 160 DIMENSION AGE(N,*), AAP(*) 161 PARAMETER ( DP5 = 0.5D0 ) 162C 163 DO 200 J = 1,N 164 JOFF = (J*J-J)/2 165 DO 100 I = 1,J 166 AAP(JOFF+I) = DP5 * (AGE(J,I) - AGE(I,J)) 167 100 CONTINUE 168 200 CONTINUE 169C 170 RETURN 171 END 172C /* Deck dgetsp */ 173 SUBROUTINE DGETSP(N,AGE,ASP) 174C 175C 8-Feb-1987 Hans Joergen Aa. Jensen 176C 177C Purpose: Transform from GE format to SP format, that is: 178C extract symmetric part of general matrix AGE 179C to symmetric, packed matrix ASP. 180C 181#include "implicit.h" 182 DIMENSION AGE(N,*), ASP(*) 183 PARAMETER ( DP5 = 0.5D0 ) 184#ifdef LUCI_DEBUG 185#include <priunit.h> 186#endif 187C 188 DO 200 J = 1,N 189 JOFF = (J*J-J)/2 190 DO 100 I = 1,J 191 ASP(JOFF+I) = DP5 * (AGE(I,J) + AGE(J,I)) 192#ifdef LUCI_DEBUG 193 write(lupri,*) 'joff,i,j ==> age(i,j), age(j,i,)',joff, 194 & i,j,AGE(I,J), AGE(J,I),ASP(JOFF+I) 195#endif 196 100 CONTINUE 197 200 CONTINUE 198#ifdef LUCI_DEBUG 199 write(lupri,*) 'end of DGETSP' 200#endif 201C 202 RETURN 203 END 204C /* Deck dgefsp */ 205 SUBROUTINE DGEFSP(N,AGE,ASP) 206C 207C 3-Nov-1989 Hans Joergen Aa. Jensen 208C 209C Purpose: Fold from GE format to SP format, that is: 210C ASP(ij) = AGE(I,J) + (1 - DELTAij) AGE(J,I) 211C 212#include "implicit.h" 213 DIMENSION AGE(N,*), ASP(*) 214C 215 DO 200 J = 1,N 216 JOFF = (J*J-J)/2 217 DO 100 I = 1,J-1 218 ASP(JOFF+I) = AGE(I,J) + AGE(J,I) 219 100 CONTINUE 220 ASP(JOFF+J) = AGE(J,J) 221 200 CONTINUE 222C 223 RETURN 224 END 225C /* Deck dgeasp */ 226 SUBROUTINE DGEASP(N,AGE,ASP) 227C 228C 4-Dec-1991 : = DGETSP but adds to SP matrix 229C 230C Purpose: Transform from GE format to SP format, that is: 231C extract symmetric part of general matrix AGE 232C to symmetric, packed matrix ASP. 233C 234#include "implicit.h" 235 DIMENSION AGE(N,*), ASP(*) 236 PARAMETER ( DP5 = 0.5D0 ) 237C 238 DO 200 J = 1,N 239 JOFF = (J*J-J)/2 240 DO 100 I = 1,J 241 ASP(JOFF+I) = ASP(JOFF+I) + DP5 * (AGE(I,J) + AGE(J,I)) 242 100 CONTINUE 243 200 CONTINUE 244C 245 RETURN 246 END 247C /* Deck dunfld */ 248 SUBROUTINE DUNFLD(N,ASP,AGE) 249C 250C 2-Dec-1991 Hans Agren 251C 252C Purpose: Unfold from SP format to GE format, that is: 253C AGE(I,J) = AGE(J,I) = ASP(ij)/(2.D0 - Delta(I,J)) 254C 255#include "implicit.h" 256 DIMENSION AGE(N,*), ASP(*) 257 PARAMETER ( DP5 = 0.5D0) 258C 259 DO 200 J = 1,N 260 JOFF = (J*J-J)/2 261 DO 100 I = 1,J-1 262 X = ASP(JOFF+I)*DP5 263 AGE(I,J) = X 264 AGE(J,I) = X 265 100 CONTINUE 266 AGE(J,J) = ASP(JOFF+J) 267 200 CONTINUE 268C 269 RETURN 270 END 271C /* Deck dsitsp */ 272 SUBROUTINE DSITSP(N,ASI,ASP) 273C 274C 8-Feb-1987 Hans Joergen Aa. Jensen 275C 276C Purpose: Transform from SI format to SP format, that is: 277C copy upper triangle of symmetric matrix ASI 278C to symmetric, packed matrix ASP. 279C 280#include "implicit.h" 281 DIMENSION ASI(N,*), ASP(*) 282C 283 DO 200 J = 1,N 284 JOFF = (J*J-J)/2 285 DO 100 I = 1,J 286 ASP(JOFF+I) = ASI(I,J) 287 100 CONTINUE 288 200 CONTINUE 289C 290 RETURN 291 END 292C /* Deck dsifsp */ 293 SUBROUTINE DSIFSP(N,ASI,ASP) 294C 295C 3-Nov-1989 Hans Joergen Aa. Jensen 296C 297C Purpose: Fold from SI format to SP format, that is: 298C ASP(ij) = ASI(I,J) + (1 - DELTAij) ASI(J,I) 299C = (2 - DELTAij) * ASI(I,J) 300C 301#include "implicit.h" 302 DIMENSION ASI(N,*), ASP(*) 303 PARAMETER (D2 = 2.0D0) 304C 305 DO 200 J = 1,N 306 JOFF = (J*J-J)/2 307 DO 100 I = 1,J-1 308 ASP(JOFF+I) = D2*ASI(I,J) 309 100 CONTINUE 310 ASP(JOFF+J) = ASI(J,J) 311 200 CONTINUE 312C 313 RETURN 314 END 315C /* Deck dgetsi */ 316 SUBROUTINE DGETSI(N,AGE,ASI) 317C 318C 3-Nov-1997 Hans Joergen Aa. Jensen 319C 320C Purpose: extract symmetric part of general matrix AGE 321C to symmetric matrix ASI. 322C Can be used in-place ( CALL DGETSI(N,A,A) ) 323C 324#include "implicit.h" 325 DIMENSION AGE(N,*), ASI(N,*) 326 PARAMETER ( DP5 = 0.5D0 ) 327C 328 DO 200 J = 1,N 329 ASI(J,J) = AGE(J,J) 330 DO 100 I = 1,J-1 331 ASI(J,I) = DP5 * (AGE(I,J) + AGE(J,I)) 332 ASI(I,J) = ASI(J,I) 333 100 CONTINUE 334 200 CONTINUE 335C 336 RETURN 337 END 338C /* Deck dgetrn */ 339 SUBROUTINE DGETRN(AGE,NROWA,NRDIMA) 340C 341C Replace AGE by AGE(transposed) 342C 343C 3-Apr-1987 HJAaJ 344C 900108-hjaaj: block with NBLK for reduced paging 345C when virtual memory 346C new name 971103-hjaaj (old name DGETRS was same as 347C a routine in LAPACK for solving linear equations; 348C when linking with complib on SGI/IRIX the LAPACK routine 349C was loaded instead of this one). 350C 351#include "implicit.h" 352 DIMENSION AGE(NRDIMA,*) 353 PARAMETER (NBLK = 128) 354 DO 400 JBLK = 1,NROWA,NBLK 355 JEND = MIN(NROWA,JBLK-1+NBLK) 356 DO 300 IBLK = 1,JBLK,NBLK 357 IEND = MIN(NROWA,IBLK-1+NBLK) 358 DO 200 J = JBLK,JEND 359 DO 100 I = IBLK,MIN(J-1,IEND) 360 SWAP = AGE(I,J) 361 AGE(I,J) = AGE(J,I) 362 AGE(J,I) = SWAP 363 100 CONTINUE 364 200 CONTINUE 365 300 CONTINUE 366 400 CONTINUE 367 RETURN 368 END 369 370#if !defined (VAR_DXML) 371C DXML library on DEC-ALPHA 372 373C Block 3: 374C DSUM, IDAMIN, IDMAX, IDMIN (supplement DASUM and IDAMAX) 375C /* Deck dsum */ 376 FUNCTION DSUM(N,DA,INCA) 377C 378C Sums elements of a vector. 379C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. 380C 30-Apr-1984 -- hjaaj -- based on DDOT from LINPACK 381C DDOT: JACK DONGARRA, LINPACK, 3/11/78. 382C 383#include "implicit.h" 384C 385 DIMENSION DA(*) 386 PARAMETER ( D0 = 0.0D0 ) 387C 388 IF (N.LE.0) THEN 389 DSUM = D0 390 RETURN 391 END IF 392 DTEMP = D0 393 IF(INCA.EQ.1)GO TO 20 394C 395C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 396C NOT EQUAL TO 1 397C 398 IA = 1 399 IF(INCA.LT.0)IA = (-N+1)*INCA + 1 400 DO 10 I = 1,N 401 DTEMP = DTEMP + DA(IA) 402 IA = IA + INCA 403 10 CONTINUE 404 DSUM = DTEMP 405 RETURN 406C 407C CODE FOR BOTH INCREMENTS EQUAL TO 1 408C 409C 410C CLEAN-UP LOOP 411C 412 20 M = MOD(N,5) 413 IF( M .EQ. 0 ) GO TO 40 414 DO 30 I = 1,M 415 DTEMP = DTEMP + DA(I) 416 30 CONTINUE 417 IF( N .LT. 5 ) GO TO 60 418 40 MP1 = M + 1 419 DO 50 I = MP1,N,5 420 DTEMP = DTEMP + DA(I) + DA(I + 1) 421 * + DA(I + 2) + DA(I + 3) + DA(I + 4) 422 50 CONTINUE 423 60 DSUM = DTEMP 424 RETURN 425 END 426C /* Deck idamin */ 427 INTEGER FUNCTION IDAMIN(N,DX,INCX) 428C 429C FINDS THE INDEX OF ELEMENT HAVING MIN. ABSOLUTE VALUE. 430C 890927-Hans Joergen Aa. Jensen 431C Based on IDAMAX by 432C JACK DONGARRA, LINPACK, 3/11/78. 433C 434#include "implicit.h" 435C 436 DIMENSION DX(*) 437C 438 IF( N .LT. 1 ) THEN 439 IDAMIN = 0 440 RETURN 441 END IF 442 IDAMIN = 1 443 IF(N.EQ.1)RETURN 444 IF(INCX.EQ.1)GO TO 20 445C 446C CODE FOR INCREMENT NOT EQUAL TO 1 447C 448 IX = 1 449 DMIN = DABS(DX(1)) 450 IX = IX + INCX 451 DO 10 I = 2,N 452 IF(DABS(DX(IX)).GE.DMIN) GO TO 5 453 IDAMIN = I 454 DMIN = DABS(DX(IX)) 455 5 IX = IX + INCX 456 10 CONTINUE 457 RETURN 458C 459C CODE FOR INCREMENT EQUAL TO 1 460C 461 20 DMIN = DABS(DX(1)) 462 DO 30 I = 2,N 463 IF(DABS(DX(I)).GE.DMIN) GO TO 30 464 IDAMIN = I 465 DMIN = DABS(DX(I)) 466 30 CONTINUE 467 RETURN 468 END 469C /* Deck idmax */ 470 INTEGER FUNCTION IDMAX(N,DX,INCX) 471C 472C FINDS THE INDEX OF ELEMENT HAVING MAX. VALUE. 473C 890105 hjaaj, based on IDAMAX by JACK DONGARRA, LINPACK, 3/11/78. 474C 475#include "implicit.h" 476C 477 DIMENSION DX(*) 478C 479 IF( N .LT. 1 ) THEN 480 IDMAX = 0 481 RETURN 482 END IF 483 IDMAX = 1 484 IF(N.EQ.1)RETURN 485 IF(INCX.EQ.1)GO TO 20 486C 487C CODE FOR INCREMENT NOT EQUAL TO 1 488C 489 IX = 1 490 DMAX = DX(1) 491 IX = IX + INCX 492 DO 10 I = 2,N 493 IF(DX(IX).LE.DMAX) GO TO 5 494 IDMAX = I 495 DMAX = DX(IX) 496 5 IX = IX + INCX 497 10 CONTINUE 498 RETURN 499C 500C CODE FOR INCREMENT EQUAL TO 1 501C 502 20 DMAX = DX(1) 503 DO 30 I = 2,N 504 IF(DX(I).LE.DMAX) GO TO 30 505 IDMAX = I 506 DMAX = DX(I) 507 30 CONTINUE 508 RETURN 509 END 510C /* Deck idmin */ 511 INTEGER FUNCTION IDMIN(N,DX,INCX) 512C 513C FINDS THE INDEX OF ELEMENT HAVING MIN. VALUE. 514C 890105 hjaaj, based on IDAMAX by JACK DONGARRA, LINPACK, 3/11/78. 515C 516#include "implicit.h" 517C 518 DIMENSION DX(*) 519C 520 IF( N .LT. 1 ) THEN 521 IDMIN = 0 522 RETURN 523 END IF 524 525 IDMIN = 1 526 IF(N.EQ.1)RETURN 527 IF(INCX.EQ.1)GO TO 20 528C 529C CODE FOR INCREMENT NOT EQUAL TO 1 530C 531 IX = 1 532 DMIN = DX(1) 533 IX = IX + INCX 534 DO 10 I = 2,N 535 IF(DX(IX).GE.DMIN) GO TO 5 536 IDMIN = I 537 DMIN = DX(IX) 538 5 IX = IX + INCX 539 10 CONTINUE 540 RETURN 541C 542C CODE FOR INCREMENT EQUAL TO 1 543C 544 20 DMIN = DX(1) 545 DO 30 I = 2,N 546 IF(DX(I).GE.DMIN) GO TO 30 547 IDMIN = I 548 DMIN = DX(I) 549 30 CONTINUE 550 RETURN 551 END 552#endif 553 554C Block 4: "iblas": ICOPY, ISCAL, ISWAP 555 556C /* Deck iblas */ 557 558 SUBROUTINE ICOPY(N,IX,INCX,IY,INCY) 559C 560C COPY integer IX TO integer IY. 561C FOR I = 0 TO N-1, COPY IX(LX+I*INCX) TO IY(LY+I*INCY), 562C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS 563C DEFINED IN A SIMILAR WAY USING INCY. 564C 565C (860516 - hjaaj - based on BLAS DCOPY) 566C 567 INTEGER IX(*),IY(*) 568 IF(N.LE.0)RETURN 569 IF(INCX.EQ.INCY) THEN 570 IF(INCX.EQ.1) GOTO 20 571 IF(INCX.GT.1) GOTO 60 572 END IF 573 5 CONTINUE 574C 575C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 576C 577 JX = 1 578 JY = 1 579 IF(INCX.LT.0)JX = (-N+1)*INCX + 1 580 IF(INCY.LT.0)JY = (-N+1)*INCY + 1 581 DO 10 I = 1,N 582 IY(JY) = IX(JX) 583 JX = JX + INCX 584 JY = JY + INCY 585 10 CONTINUE 586 RETURN 587C 588C CODE FOR BOTH INCREMENTS EQUAL TO 1 589C 590C 591C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. 592C 593 20 M = MOD(N,7) 594 IF( M .EQ. 0 ) GO TO 40 595 DO 30 I = 1,M 596 IY(I) = IX(I) 597 30 CONTINUE 598 IF( N .LT. 7 ) RETURN 599 40 MP1 = M + 1 600 DO 50 I = MP1,N,7 601 IY(I) = IX(I) 602 IY(I + 1) = IX(I + 1) 603 IY(I + 2) = IX(I + 2) 604 IY(I + 3) = IX(I + 3) 605 IY(I + 4) = IX(I + 4) 606 IY(I + 5) = IX(I + 5) 607 IY(I + 6) = IX(I + 6) 608 50 CONTINUE 609 RETURN 610C 611C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. 612C 613 60 CONTINUE 614 NS=N*INCX 615 DO 70 I=1,NS,INCX 616 IY(I) = IX(I) 617 70 CONTINUE 618 RETURN 619 END 620 SUBROUTINE ISCAL(N,IA,IX,INCX) 621C 622C Scale integer vector IX with IA 623C FOR I = 0 TO N-1, SCALE IX(LX+I*INCX) WITH IA 624C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N 625C 626C (901219 - hjaaj - based on ICOPY) 627CC 5) DSPSOL, DGESOL, DGEINV : driver routines calling LINPACK 628 INTEGER IX(*) 629 IF(N.LE.0)RETURN 630 IF(INCX.EQ.1) GOTO 20 631 IF(INCX.GT.1) GOTO 60 632 5 CONTINUE 633C 634C CODE FOR NONPOSITIVE INCREMENT. 635C 636 JX = (-N+1)*INCX + 1 637 DO 10 I = 1,N 638 IX(JX) = IA*IX(JX) 639 JX = JX + INCX 640 10 CONTINUE 641 RETURN 642C 643C CODE FOR INCREMENT EQUAL TO 1 644C 645C 646C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. 647C 648 20 M = MOD(N,7) 649 IF( M .EQ. 0 ) GO TO 40 650 DO 30 I = 1,M 651 IX(I) = IA*IX(I) 652 30 CONTINUE 653 IF( N .LT. 7 ) RETURN 654 40 MP1 = M + 1 655 DO 50 I = MP1,N,7 656 IX(I) = IA*IX(I) 657 IX(I + 1) = IA*IX(I + 1) 658 IX(I + 2) = IA*IX(I + 2) 659 IX(I + 3) = IA*IX(I + 3) 660 IX(I + 4) = IA*IX(I + 4) 661 IX(I + 5) = IA*IX(I + 5) 662 IX(I + 6) = IA*IX(I + 6) 663 50 CONTINUE 664 RETURN 665C 666C CODE FOR POSITIVE, NONUNIT INCREMENT. 667C 668 60 CONTINUE 669 NS=N*INCX 670 DO 70 I=1,NS,INCX 671 IX(I) = IA*IX(I) 672 70 CONTINUE 673 RETURN 674 END 675 SUBROUTINE ISWAP(N,IX,INCX,IY,INCY) 676C 677C Swap integer arrays IX and IY. 678C FOR I = 0 TO N-1, SWAP IX(LX+I*INCX) WITH IY(LY+I*INCY), 679C WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS 680C DEFINED IN A SIMILAR WAY USING INCY. 681C 682C (901219 - hjaaj - based on ICOPY) 683C 684 INTEGER IX(*),IY(*) 685 IF(N.LE.0)RETURN 686 IF(INCX.EQ.INCY .AND. INCX .GT. 0) GO TO 60 687C 688C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 689C 690 JX = 1 691 JY = 1 692 IF(INCX.LT.0)JX = (-N+1)*INCX + 1 693 IF(INCY.LT.0)JY = (-N+1)*INCY + 1 694 DO 10 I = 1,N 695 IHOLD = IY(JY) 696 IY(JY) = IX(JX) 697 IX(JY) = IHOLD 698 JX = JX + INCX 699 JY = JY + INCY 700 10 CONTINUE 701 RETURN 702C 703C CODE FOR EQUAL, POSITIVE INCREMENTS. 704C 705 60 CONTINUE 706 NS=N*INCX 707 DO 70 I=1,NS,INCX 708 IHOLD = IY(I) 709 IY(I) = IX(I) 710 IX(I) = IHOLD 711 70 CONTINUE 712 RETURN 713 END 714 715C Block 5: other routines 716C DGEZERO, DSPZERO, NDXGTA, DUNIT, DZERO, ISUM, IZERO 717C 718 719 SUBROUTINE DGEZERO(N,AGE,NRSTA,NREND,NCSTA,NCEND) 720C 721C Oct. 09, Hans Joergen Aa. Jensen 722C 723C Zero block AGE(NRSTA:NREND,NCSTA:NCEND) 724C 725 IMPLICIT NONE 726 INTEGER N, NRSTA, NREND, NCSTA, NCEND 727 REAL*8 AGE(N,N) 728 INTEGER IC, IR 729 730 DO IC = NCSTA, NCEND 731 DO IR = NRSTA, NREND 732 AGE(IR,IC) = 0.0D0 733 END DO 734 END DO 735 RETURN 736 END 737 738 SUBROUTINE DSPZERO(N,ASP,NRSTA,NREND,NCSTA,NCEND) 739C 740C Oct. 09, Hans Joergen Aa. Jensen 741C 742C Zero block ASP(NRSTA:NREND,NCSTA:NCEND), 743C however, ASP is Symmetric Packed (triangular packed) 744C 745 IMPLICIT NONE 746 REAL*8 ASP(*) 747 INTEGER N, NRSTA, NREND, NCSTA, NCEND 748 INTEGER IC, IR, IROFF, MCEND, ICOFF, MREND 749 750C First the part placed in lower triangle of the full matrix 751 752 753 IF (NCSTA .LE. NREND) THEN 754 DO IR = NRSTA, NREND 755 IROFF = (IR*IR-IR)/2 756 MCEND = MIN(NCEND, IR) 757 DO IC = NCSTA, MCEND 758 ASP(IROFF+IC) = 0.0D0 759 END DO 760 END DO 761 END IF 762 763C and then the part placed in uppper triangle 764 765 IF (NRSTA .LE. NCEND) THEN 766 DO IC = NCSTA, NCEND 767 ICOFF = (IC*IC-IC)/2 768 MREND = MIN(NREND, IC) 769 DO IR = NRSTA, MREND 770 ASP(ICOFF+IR) = 0.0D0 771 END DO 772 END DO 773 END IF 774 775 RETURN 776 END 777 778C /* Deck ndxgta */ 779 INTEGER FUNCTION NDXGTA(N,A,DX,INCX) 780C 781C 900319-hjaaj 782C 783C Return number of elements with absolute value .gt. A 784C 785#include "implicit.h" 786 DIMENSION DX(N) 787 IF (A .LT. 0.0D0) THEN 788 NUM = N 789 ELSE IF (INCX .EQ. 1) THEN 790 NUM = 0 791 DO 200 I = 1,N 792 IF (ABS(DX(I)) .GT. A) NUM = NUM + 1 793 200 CONTINUE 794 ELSE 795 NUM = 0 796 IF (INCX.GT.0) THEN 797 IX = 1 - INCX 798 ELSE 799 IX = 1 - N*INCX 800 END IF 801 DO 300 I = 1,N 802 IF (ABS(DX(IX+I*INCX)) .GT. A) NUM = NUM + 1 803 300 CONTINUE 804 END IF 805 NDXGTA = NUM 806 RETURN 807 END 808C /* Deck dunit */ 809 SUBROUTINE DUNIT(A,N) 810C 811C SUBROUTINE DUNIT SETS THE REAL SQUARE MATRIX A EQUAL 812C TO A UNIT MATRIX. 813C /VER 2/ 14-Sep-1983 hjaaj 814C 815#include "implicit.h" 816 DIMENSION A(*) 817 PARAMETER (D1=1.0D00, D0=0.0D00) 818C 819 NN = N*N 820 DO 100 I = 1,NN 821 A(I) = D0 822 100 CONTINUE 823 N1 = N + 1 824 DO 200 I = 1,NN,N1 825 A(I) = D1 826 200 CONTINUE 827 RETURN 828 END 829C /* Deck dzero */ 830 SUBROUTINE DZERO(DX,LENGTH) 831#include "implicit.h" 832C 833C Last revision 5-May-1984 by Hans Jorgen Aa. Jensen 834C 835C Subroutine DZERO sets a real array of length *LENGTH* 836C to zero. 837C................................................................... 838 DIMENSION DX(*) 839C 840 IF (LENGTH.LE.0) RETURN 841C 842 DO I = 1,LENGTH 843 DX(I) = 0.0D00 844 END DO 845C 846 RETURN 847 END 848C /* Deck isum */ 849 FUNCTION ISUM(N,IA,INCA) 850C 851C 8-Feb-1987 hjaaj 852C Sums elements of a integer vector. 853C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. 854C -- based on DDOT from LINPACK 855C DDOT: JACK DONGARRA, LINPACK, 3/11/78. 856C 857 INTEGER ISUM, IA(*), ITEMP 858 INTEGER I,INCA,JA,M,MP1,N 859C 860 IF(N.LE.0) THEN 861 ISUM = 0 862 RETURN 863 END IF 864 865 ITEMP = 0 866 IF(INCA.EQ.1)GO TO 20 867C 868C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 869C NOT EQUAL TO 1 870C 871 JA = 1 872 IF(INCA.LT.0)JA = (-N+1)*INCA + 1 873 DO 10 I = 1,N 874 ITEMP = ITEMP + IA(JA) 875 JA = JA + INCA 876 10 CONTINUE 877 ISUM = ITEMP 878 RETURN 879C 880C CODE FOR BOTH INCREMENTS EQUAL TO 1 881C 882C 883C CLEAN-UP LOOP 884C 885 20 M = MOD(N,5) 886 IF( M .EQ. 0 ) GO TO 40 887 DO 30 I = 1,M 888 ITEMP = ITEMP + IA(I) 889 30 CONTINUE 890 IF( N .LT. 5 ) GO TO 60 891 40 MP1 = M + 1 892 DO 50 I = MP1,N,5 893 ITEMP = ITEMP + IA(I) + IA(I + 1) 894 * + IA(I + 2) + IA(I + 3) + IA(I + 4) 895 50 CONTINUE 896 60 ISUM = ITEMP 897 RETURN 898 END 899C /* Deck izero */ 900 SUBROUTINE IZERO(INTVEC,LENGTH) 901C................................................................... 902C Written 5-May-1984 by Hans Jorgen Aa. Jensen 903C 904C Subroutine IZERO sets an integer array of length *LENGTH* 905C to zero. 906C................................................................... 907 INTEGER LENGTH, INTVEC(*) 908C 909 IF (LENGTH.LE.0) RETURN 910C 911 DO I=1,LENGTH 912 INTVEC(I) = 0 913 END DO 914C 915 END 916 917 918C Block 6: 919C DGEINV, DGESOL, DSPSOL, DSPSLI : driver routines calling LINPACK 920 921C /* Deck dgeinv */ 922 SUBROUTINE DGEINV(N,A,AINV,IPVT,WRK,INFO) 923C 924C 850618-HJAAJ 925C 926C Call Linpack routines to calculate the inverse of a 927C general matrix A. 928C 929 INTEGER N, IPVT(*), N2 930 REAL*8 A(*), AINV(*), WRK(*), DET(2) 931C 932 N2 = N*N 933 CALL DCOPY(N2,A,1,AINV,1) 934 CALL DGEFA(AINV,N,N,IPVT,INFO) 935 IF (INFO .EQ. 0) CALL DGEDI(AINV,N,N,IPVT,DET,WRK,01) 936 RETURN 937 END 938C /* Deck dgesol */ 939 SUBROUTINE DGESOL (NSIM,N,LDA,LDB,A,B,KPVT,INFO) 940C 941C Written 22-Mar-1985 Hans Joergen Aa. Jensen 942C No revisions. 943C 944C Purpose: 945C Solve the NSIM simultaneous eqautions: 946C 947C B(n,nsim) := A(n,n) inverse * B(n,nsim) 948C 949C using LINPACK routines DGEFA and DGESL. 950C 951C A is the matrix (general, non-singular) 952C KPVT is a scratch array of length N. 953C 954#include "implicit.h" 955 DIMENSION A(LDA,*),B(LDB,*),KPVT(*) 956C 957 CALL DGEFA (A,LDA,N,KPVT,INFO) 958 IF (INFO.NE.0) RETURN 959C 960 DO 100 J = 1,NSIM 961 CALL DGESL (A,LDA,N,KPVT,B(1,J),0) 962 100 CONTINUE 963C 964 RETURN 965 END 966C /* Deck dspsol */ 967 SUBROUTINE DSPSOL (N,NSIM,AP,B,KPVT,INFO) 968C 969C Written 8-Feb-1985 Hans Joergen Aa. Jensen 970C No revisions. 971C 972C Purpose: 973C Solve the NSIM simultaneous eqautions: 974C 975C B(n,nsim) := A(n,n) inverse * B(n,nsim) 976C 977C AP is A in packed format. 978C KPVT is a scratch array of length N. 979C 980#include "implicit.h" 981 DIMENSION AP(*),B(N,*),KPVT(*) 982C 983 CALL DSPFA (AP,N,KPVT,INFO) 984 IF (INFO.NE.0) RETURN 985C 986 DO 100 J = 1,NSIM 987 CALL DSPSL (AP,N,KPVT,B(1,J)) 988 100 CONTINUE 989C 990 RETURN 991 END 992C /* Deck dspsli */ 993 SUBROUTINE DSPSLI (N,NSIM,AP,B,KPVT,INFO,DET,INERT) 994C 995C Written 24-Feb-1989 Hans Joergen Aa. Jensen, based on DSPSOL. 996C No revisions. 997C 998C Purpose: 999C Solve the NSIM simultaneous eqautions: 1000C 1001C B(n,nsim) := A(n,n) inverse * B(n,nsim) 1002C 1003C AP is A in packed format. 1004C KPVT is a scratch array of length N. 1005C 1006#include "implicit.h" 1007 DIMENSION AP(*),B(N,*),KPVT(*),DET(2),INERT(3) 1008C 1009 CALL DSPFA (AP,N,KPVT,INFO) 1010 IF (INFO.NE.0) RETURN 1011C 1012 CALL DSPDI(AP,N,KPVT,DET,INERT,DUMMY,110) 1013C CALL DSPDI(AP,N,KPVT,DET,INERT,WORK,JOB) 1014C 1015 DO 100 J = 1,NSIM 1016 CALL DSPSL (AP,N,KPVT,B(1,J)) 1017 100 CONTINUE 1018C 1019 RETURN 1020 END 1021C Block 7: BLAS and LAPACK routines missing for ESSL 1022C LSAME; ILAENV, DSPTRI, DSPCON 1023#if defined (VAR_ESSL) 1024 LOGICAL FUNCTION LSAME(CA,CB) 1025* 1026* -- Reference BLAS level1 routine (version 3.1) -- 1027* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- 1028* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 1029* November 2011 1030* 1031* .. Scalar Arguments .. 1032 CHARACTER CA,CB 1033* .. 1034* 1035* ===================================================================== 1036* 1037* .. Intrinsic Functions .. 1038 INTRINSIC ICHAR 1039* .. 1040* .. Local Scalars .. 1041 INTEGER INTA,INTB,ZCODE 1042* .. 1043* 1044* Test if the characters are equal 1045* 1046 LSAME = CA .EQ. CB 1047 IF (LSAME) RETURN 1048* 1049* Now test for equivalence if both characters are alphabetic. 1050* 1051 ZCODE = ICHAR('Z') 1052* 1053* Use 'Z' rather than 'A' so that ASCII can be detected on Prime 1054* machines, on which ICHAR returns a value with bit 8 set. 1055* ICHAR('A') on Prime machines returns 193 which is the same as 1056* ICHAR('A') on an EBCDIC machine. 1057* 1058 INTA = ICHAR(CA) 1059 INTB = ICHAR(CB) 1060* 1061 IF (ZCODE.EQ.90 .OR. ZCODE.EQ.122) THEN 1062* 1063* ASCII is assumed - ZCODE is the ASCII code of either lower or 1064* upper case 'Z'. 1065* 1066 IF (INTA.GE.97 .AND. INTA.LE.122) INTA = INTA - 32 1067 IF (INTB.GE.97 .AND. INTB.LE.122) INTB = INTB - 32 1068* 1069 ELSE IF (ZCODE.EQ.233 .OR. ZCODE.EQ.169) THEN 1070* 1071* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or 1072* upper case 'Z'. 1073* 1074 IF (INTA.GE.129 .AND. INTA.LE.137 .OR. 1075 + INTA.GE.145 .AND. INTA.LE.153 .OR. 1076 + INTA.GE.162 .AND. INTA.LE.169) INTA = INTA + 64 1077 IF (INTB.GE.129 .AND. INTB.LE.137 .OR. 1078 + INTB.GE.145 .AND. INTB.LE.153 .OR. 1079 + INTB.GE.162 .AND. INTB.LE.169) INTB = INTB + 64 1080* 1081 ELSE IF (ZCODE.EQ.218 .OR. ZCODE.EQ.250) THEN 1082* 1083* ASCII is assumed, on Prime machines - ZCODE is the ASCII code 1084* plus 128 of either lower or upper case 'Z'. 1085* 1086 IF (INTA.GE.225 .AND. INTA.LE.250) INTA = INTA - 32 1087 IF (INTB.GE.225 .AND. INTB.LE.250) INTB = INTB - 32 1088 END IF 1089 LSAME = INTA .EQ. INTB 1090* 1091* RETURN 1092* 1093* End of LSAME 1094* 1095 END 1096 INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 ) 1097* 1098* -- LAPACK auxiliary routine (version 3.4.0) -- 1099* -- LAPACK is a software package provided by Univ. of Tennessee, -- 1100* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 1101* November 2011 1102* 1103* .. Scalar Arguments .. 1104 CHARACTER*( * ) NAME, OPTS 1105 INTEGER ISPEC, N1, N2, N3, N4 1106* .. 1107* 1108* ===================================================================== 1109* 1110* .. Local Scalars .. 1111 INTEGER I, IC, IZ, NB, NBMIN, NX 1112 LOGICAL CNAME, SNAME 1113 CHARACTER C1*1, C2*2, C4*2, C3*3, SUBNAM*6 1114* .. 1115* .. Intrinsic Functions .. 1116 INTRINSIC CHAR, ICHAR, INT, MIN, REAL 1117* .. 1118* .. External Functions .. 1119 INTEGER IEEECK, IPARMQ 1120 EXTERNAL IEEECK, IPARMQ 1121* .. 1122* .. Executable Statements .. 1123* 1124 GO TO ( 10, 10, 10, 80, 90, 100, 110, 120, 1125 $ 130, 140, 150, 160, 160, 160, 160, 160 )ISPEC 1126* 1127* Invalid value for ISPEC 1128* 1129 ILAENV = -1 1130 RETURN 1131* 1132 10 CONTINUE 1133* 1134* Convert NAME to upper case if the first character is lower case. 1135* 1136 ILAENV = 1 1137 SUBNAM = NAME 1138 IC = ICHAR( SUBNAM( 1: 1 ) ) 1139 IZ = ICHAR( 'Z' ) 1140 IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN 1141* 1142* ASCII character set 1143* 1144 IF( IC.GE.97 .AND. IC.LE.122 ) THEN 1145 SUBNAM( 1: 1 ) = CHAR( IC-32 ) 1146 DO 20 I = 2, 6 1147 IC = ICHAR( SUBNAM( I: I ) ) 1148 IF( IC.GE.97 .AND. IC.LE.122 ) 1149 $ SUBNAM( I: I ) = CHAR( IC-32 ) 1150 20 CONTINUE 1151 END IF 1152* 1153 ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN 1154* 1155* EBCDIC character set 1156* 1157 IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. 1158 $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. 1159 $ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN 1160 SUBNAM( 1: 1 ) = CHAR( IC+64 ) 1161 DO 30 I = 2, 6 1162 IC = ICHAR( SUBNAM( I: I ) ) 1163 IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. 1164 $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. 1165 $ ( IC.GE.162 .AND. IC.LE.169 ) )SUBNAM( I: 1166 $ I ) = CHAR( IC+64 ) 1167 30 CONTINUE 1168 END IF 1169* 1170 ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN 1171* 1172* Prime machines: ASCII+128 1173* 1174 IF( IC.GE.225 .AND. IC.LE.250 ) THEN 1175 SUBNAM( 1: 1 ) = CHAR( IC-32 ) 1176 DO 40 I = 2, 6 1177 IC = ICHAR( SUBNAM( I: I ) ) 1178 IF( IC.GE.225 .AND. IC.LE.250 ) 1179 $ SUBNAM( I: I ) = CHAR( IC-32 ) 1180 40 CONTINUE 1181 END IF 1182 END IF 1183* 1184 C1 = SUBNAM( 1: 1 ) 1185 SNAME = C1.EQ.'S' .OR. C1.EQ.'D' 1186 CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' 1187 IF( .NOT.( CNAME .OR. SNAME ) ) 1188 $ RETURN 1189 C2 = SUBNAM( 2: 3 ) 1190 C3 = SUBNAM( 4: 6 ) 1191 C4 = C3( 2: 3 ) 1192* 1193 GO TO ( 50, 60, 70 )ISPEC 1194* 1195 50 CONTINUE 1196* 1197* ISPEC = 1: block size 1198* 1199* In these examples, separate code is provided for setting NB for 1200* real and complex. We assume that NB will take the same value in 1201* single or double precision. 1202* 1203 NB = 1 1204* 1205 IF( C2.EQ.'GE' ) THEN 1206 IF( C3.EQ.'TRF' ) THEN 1207 IF( SNAME ) THEN 1208 NB = 64 1209 ELSE 1210 NB = 64 1211 END IF 1212 ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. 1213 $ C3.EQ.'QLF' ) THEN 1214 IF( SNAME ) THEN 1215 NB = 32 1216 ELSE 1217 NB = 32 1218 END IF 1219 ELSE IF( C3.EQ.'HRD' ) THEN 1220 IF( SNAME ) THEN 1221 NB = 32 1222 ELSE 1223 NB = 32 1224 END IF 1225 ELSE IF( C3.EQ.'BRD' ) THEN 1226 IF( SNAME ) THEN 1227 NB = 32 1228 ELSE 1229 NB = 32 1230 END IF 1231 ELSE IF( C3.EQ.'TRI' ) THEN 1232 IF( SNAME ) THEN 1233 NB = 64 1234 ELSE 1235 NB = 64 1236 END IF 1237 END IF 1238 ELSE IF( C2.EQ.'PO' ) THEN 1239 IF( C3.EQ.'TRF' ) THEN 1240 IF( SNAME ) THEN 1241 NB = 64 1242 ELSE 1243 NB = 64 1244 END IF 1245 END IF 1246 ELSE IF( C2.EQ.'SY' ) THEN 1247 IF( C3.EQ.'TRF' ) THEN 1248 IF( SNAME ) THEN 1249 NB = 64 1250 ELSE 1251 NB = 64 1252 END IF 1253 ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN 1254 NB = 32 1255 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN 1256 NB = 64 1257 END IF 1258 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 1259 IF( C3.EQ.'TRF' ) THEN 1260 NB = 64 1261 ELSE IF( C3.EQ.'TRD' ) THEN 1262 NB = 32 1263 ELSE IF( C3.EQ.'GST' ) THEN 1264 NB = 64 1265 END IF 1266 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 1267 IF( C3( 1: 1 ).EQ.'G' ) THEN 1268 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1269 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1270 $ THEN 1271 NB = 32 1272 END IF 1273 ELSE IF( C3( 1: 1 ).EQ.'M' ) THEN 1274 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1275 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1276 $ THEN 1277 NB = 32 1278 END IF 1279 END IF 1280 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 1281 IF( C3( 1: 1 ).EQ.'G' ) THEN 1282 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1283 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1284 $ THEN 1285 NB = 32 1286 END IF 1287 ELSE IF( C3( 1: 1 ).EQ.'M' ) THEN 1288 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1289 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1290 $ THEN 1291 NB = 32 1292 END IF 1293 END IF 1294 ELSE IF( C2.EQ.'GB' ) THEN 1295 IF( C3.EQ.'TRF' ) THEN 1296 IF( SNAME ) THEN 1297 IF( N4.LE.64 ) THEN 1298 NB = 1 1299 ELSE 1300 NB = 32 1301 END IF 1302 ELSE 1303 IF( N4.LE.64 ) THEN 1304 NB = 1 1305 ELSE 1306 NB = 32 1307 END IF 1308 END IF 1309 END IF 1310 ELSE IF( C2.EQ.'PB' ) THEN 1311 IF( C3.EQ.'TRF' ) THEN 1312 IF( SNAME ) THEN 1313 IF( N2.LE.64 ) THEN 1314 NB = 1 1315 ELSE 1316 NB = 32 1317 END IF 1318 ELSE 1319 IF( N2.LE.64 ) THEN 1320 NB = 1 1321 ELSE 1322 NB = 32 1323 END IF 1324 END IF 1325 END IF 1326 ELSE IF( C2.EQ.'TR' ) THEN 1327 IF( C3.EQ.'TRI' ) THEN 1328 IF( SNAME ) THEN 1329 NB = 64 1330 ELSE 1331 NB = 64 1332 END IF 1333 END IF 1334 ELSE IF( C2.EQ.'LA' ) THEN 1335 IF( C3.EQ.'UUM' ) THEN 1336 IF( SNAME ) THEN 1337 NB = 64 1338 ELSE 1339 NB = 64 1340 END IF 1341 END IF 1342 ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN 1343 IF( C3.EQ.'EBZ' ) THEN 1344 NB = 1 1345 END IF 1346 END IF 1347 ILAENV = NB 1348 RETURN 1349* 1350 60 CONTINUE 1351* 1352* ISPEC = 2: minimum block size 1353* 1354 NBMIN = 2 1355 IF( C2.EQ.'GE' ) THEN 1356 IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. C3.EQ. 1357 $ 'QLF' ) THEN 1358 IF( SNAME ) THEN 1359 NBMIN = 2 1360 ELSE 1361 NBMIN = 2 1362 END IF 1363 ELSE IF( C3.EQ.'HRD' ) THEN 1364 IF( SNAME ) THEN 1365 NBMIN = 2 1366 ELSE 1367 NBMIN = 2 1368 END IF 1369 ELSE IF( C3.EQ.'BRD' ) THEN 1370 IF( SNAME ) THEN 1371 NBMIN = 2 1372 ELSE 1373 NBMIN = 2 1374 END IF 1375 ELSE IF( C3.EQ.'TRI' ) THEN 1376 IF( SNAME ) THEN 1377 NBMIN = 2 1378 ELSE 1379 NBMIN = 2 1380 END IF 1381 END IF 1382 ELSE IF( C2.EQ.'SY' ) THEN 1383 IF( C3.EQ.'TRF' ) THEN 1384 IF( SNAME ) THEN 1385 NBMIN = 8 1386 ELSE 1387 NBMIN = 8 1388 END IF 1389 ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN 1390 NBMIN = 2 1391 END IF 1392 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 1393 IF( C3.EQ.'TRD' ) THEN 1394 NBMIN = 2 1395 END IF 1396 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 1397 IF( C3( 1: 1 ).EQ.'G' ) THEN 1398 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1399 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1400 $ THEN 1401 NBMIN = 2 1402 END IF 1403 ELSE IF( C3( 1: 1 ).EQ.'M' ) THEN 1404 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1405 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1406 $ THEN 1407 NBMIN = 2 1408 END IF 1409 END IF 1410 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 1411 IF( C3( 1: 1 ).EQ.'G' ) THEN 1412 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1413 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1414 $ THEN 1415 NBMIN = 2 1416 END IF 1417 ELSE IF( C3( 1: 1 ).EQ.'M' ) THEN 1418 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1419 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1420 $ THEN 1421 NBMIN = 2 1422 END IF 1423 END IF 1424 END IF 1425 ILAENV = NBMIN 1426 RETURN 1427* 1428 70 CONTINUE 1429* 1430* ISPEC = 3: crossover point 1431* 1432 NX = 0 1433 IF( C2.EQ.'GE' ) THEN 1434 IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. C3.EQ. 1435 $ 'QLF' ) THEN 1436 IF( SNAME ) THEN 1437 NX = 128 1438 ELSE 1439 NX = 128 1440 END IF 1441 ELSE IF( C3.EQ.'HRD' ) THEN 1442 IF( SNAME ) THEN 1443 NX = 128 1444 ELSE 1445 NX = 128 1446 END IF 1447 ELSE IF( C3.EQ.'BRD' ) THEN 1448 IF( SNAME ) THEN 1449 NX = 128 1450 ELSE 1451 NX = 128 1452 END IF 1453 END IF 1454 ELSE IF( C2.EQ.'SY' ) THEN 1455 IF( SNAME .AND. C3.EQ.'TRD' ) THEN 1456 NX = 32 1457 END IF 1458 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 1459 IF( C3.EQ.'TRD' ) THEN 1460 NX = 32 1461 END IF 1462 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 1463 IF( C3( 1: 1 ).EQ.'G' ) THEN 1464 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1465 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1466 $ THEN 1467 NX = 128 1468 END IF 1469 END IF 1470 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 1471 IF( C3( 1: 1 ).EQ.'G' ) THEN 1472 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. C4.EQ. 1473 $ 'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. C4.EQ.'BR' ) 1474 $ THEN 1475 NX = 128 1476 END IF 1477 END IF 1478 END IF 1479 ILAENV = NX 1480 RETURN 1481* 1482 80 CONTINUE 1483* 1484* ISPEC = 4: number of shifts (used by xHSEQR) 1485* 1486 ILAENV = 6 1487 RETURN 1488* 1489 90 CONTINUE 1490* 1491* ISPEC = 5: minimum column dimension (not used) 1492* 1493 ILAENV = 2 1494 RETURN 1495* 1496 100 CONTINUE 1497* 1498* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) 1499* 1500 ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) 1501 RETURN 1502* 1503 110 CONTINUE 1504* 1505* ISPEC = 7: number of processors (not used) 1506* 1507 ILAENV = 1 1508 RETURN 1509* 1510 120 CONTINUE 1511* 1512* ISPEC = 8: crossover point for multishift (used by xHSEQR) 1513* 1514 ILAENV = 50 1515 RETURN 1516* 1517 130 CONTINUE 1518* 1519* ISPEC = 9: maximum size of the subproblems at the bottom of the 1520* computation tree in the divide-and-conquer algorithm 1521* (used by xGELSD and xGESDD) 1522* 1523 ILAENV = 25 1524 RETURN 1525* 1526 140 CONTINUE 1527* 1528* ISPEC = 10: ieee NaN arithmetic can be trusted not to trap 1529* 1530* ILAENV = 0 1531 ILAENV = 1 1532 IF( ILAENV.EQ.1 ) THEN 1533 ILAENV = IEEECK( 1, 0.0, 1.0 ) 1534 END IF 1535 RETURN 1536* 1537 150 CONTINUE 1538* 1539* ISPEC = 11: infinity arithmetic can be trusted not to trap 1540* 1541* ILAENV = 0 1542 ILAENV = 1 1543 IF( ILAENV.EQ.1 ) THEN 1544 ILAENV = IEEECK( 0, 0.0, 1.0 ) 1545 END IF 1546 RETURN 1547* 1548 160 CONTINUE 1549* 1550* 12 <= ISPEC <= 16: xHSEQR or one of its subroutines. 1551* 1552 ILAENV = IPARMQ( ISPEC, NAME, OPTS, N1, N2, N3, N4 ) 1553 RETURN 1554* 1555* End of ILAENV 1556* 1557 END 1558 SUBROUTINE DSPTRI( UPLO, N, AP, IPIV, WORK, INFO ) 1559* 1560* -- LAPACK computational routine (version 3.4.0) -- 1561* -- LAPACK is a software package provided by Univ. of Tennessee, -- 1562* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 1563* November 2011 1564* 1565* .. Scalar Arguments .. 1566 CHARACTER UPLO 1567 INTEGER INFO, N 1568* .. 1569* .. Array Arguments .. 1570 INTEGER IPIV( * ) 1571 DOUBLE PRECISION AP( * ), WORK( * ) 1572* .. 1573* 1574* ===================================================================== 1575* 1576* .. Parameters .. 1577 DOUBLE PRECISION ONE, ZERO 1578 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 1579* .. 1580* .. Local Scalars .. 1581 LOGICAL UPPER 1582 INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP 1583 DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP 1584* .. 1585* .. External Functions .. 1586 LOGICAL LSAME 1587 DOUBLE PRECISION DDOT 1588 EXTERNAL LSAME, DDOT 1589* .. 1590* .. External Subroutines .. 1591 EXTERNAL DCOPY, DSPMV, DSWAP, XERBLA 1592* .. 1593* .. Intrinsic Functions .. 1594 INTRINSIC ABS 1595* .. 1596* .. Executable Statements .. 1597* 1598* Test the input parameters. 1599* 1600 INFO = 0 1601 UPPER = LSAME( UPLO, 'U' ) 1602 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 1603 INFO = -1 1604 ELSE IF( N.LT.0 ) THEN 1605 INFO = -2 1606 END IF 1607 IF( INFO.NE.0 ) THEN 1608 CALL XERBLA( 'DSPTRI', -INFO ) 1609 RETURN 1610 END IF 1611* 1612* Quick return if possible 1613* 1614 IF( N.EQ.0 ) 1615 $ RETURN 1616* 1617* Check that the diagonal matrix D is nonsingular. 1618* 1619 IF( UPPER ) THEN 1620* 1621* Upper triangular storage: examine D from bottom to top 1622* 1623 KP = N*( N+1 ) / 2 1624 DO 10 INFO = N, 1, -1 1625 IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO ) 1626 $ RETURN 1627 KP = KP - INFO 1628 10 CONTINUE 1629 ELSE 1630* 1631* Lower triangular storage: examine D from top to bottom. 1632* 1633 KP = 1 1634 DO 20 INFO = 1, N 1635 IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO ) 1636 $ RETURN 1637 KP = KP + N - INFO + 1 1638 20 CONTINUE 1639 END IF 1640 INFO = 0 1641* 1642 IF( UPPER ) THEN 1643* 1644* Compute inv(A) from the factorization A = U*D*U**T. 1645* 1646* K is the main loop index, increasing from 1 to N in steps of 1647* 1 or 2, depending on the size of the diagonal blocks. 1648* 1649 K = 1 1650 KC = 1 1651 30 CONTINUE 1652* 1653* If K > N, exit from loop. 1654* 1655 IF( K.GT.N ) 1656 $ GO TO 50 1657* 1658 KCNEXT = KC + K 1659 IF( IPIV( K ).GT.0 ) THEN 1660* 1661* 1 x 1 diagonal block 1662* 1663* Invert the diagonal block. 1664* 1665 AP( KC+K-1 ) = ONE / AP( KC+K-1 ) 1666* 1667* Compute column K of the inverse. 1668* 1669 IF( K.GT.1 ) THEN 1670 CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 ) 1671 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ), 1672 $ 1 ) 1673 AP( KC+K-1 ) = AP( KC+K-1 ) - 1674 $ DDOT( K-1, WORK, 1, AP( KC ), 1 ) 1675 END IF 1676 KSTEP = 1 1677 ELSE 1678* 1679* 2 x 2 diagonal block 1680* 1681* Invert the diagonal block. 1682* 1683 T = ABS( AP( KCNEXT+K-1 ) ) 1684 AK = AP( KC+K-1 ) / T 1685 AKP1 = AP( KCNEXT+K ) / T 1686 AKKP1 = AP( KCNEXT+K-1 ) / T 1687 D = T*( AK*AKP1-ONE ) 1688 AP( KC+K-1 ) = AKP1 / D 1689 AP( KCNEXT+K ) = AK / D 1690 AP( KCNEXT+K-1 ) = -AKKP1 / D 1691* 1692* Compute columns K and K+1 of the inverse. 1693* 1694 IF( K.GT.1 ) THEN 1695 CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 ) 1696 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ), 1697 $ 1 ) 1698 AP( KC+K-1 ) = AP( KC+K-1 ) - 1699 $ DDOT( K-1, WORK, 1, AP( KC ), 1 ) 1700 AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) - 1701 $ DDOT( K-1, AP( KC ), 1, AP( KCNEXT ), 1702 $ 1 ) 1703 CALL DCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 ) 1704 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, 1705 $ AP( KCNEXT ), 1 ) 1706 AP( KCNEXT+K ) = AP( KCNEXT+K ) - 1707 $ DDOT( K-1, WORK, 1, AP( KCNEXT ), 1 ) 1708 END IF 1709 KSTEP = 2 1710 KCNEXT = KCNEXT + K + 1 1711 END IF 1712* 1713 KP = ABS( IPIV( K ) ) 1714 IF( KP.NE.K ) THEN 1715* 1716* Interchange rows and columns K and KP in the leading 1717* submatrix A(1:k+1,1:k+1) 1718* 1719 KPC = ( KP-1 )*KP / 2 + 1 1720 CALL DSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 ) 1721 KX = KPC + KP - 1 1722 DO 40 J = KP + 1, K - 1 1723 KX = KX + J - 1 1724 TEMP = AP( KC+J-1 ) 1725 AP( KC+J-1 ) = AP( KX ) 1726 AP( KX ) = TEMP 1727 40 CONTINUE 1728 TEMP = AP( KC+K-1 ) 1729 AP( KC+K-1 ) = AP( KPC+KP-1 ) 1730 AP( KPC+KP-1 ) = TEMP 1731 IF( KSTEP.EQ.2 ) THEN 1732 TEMP = AP( KC+K+K-1 ) 1733 AP( KC+K+K-1 ) = AP( KC+K+KP-1 ) 1734 AP( KC+K+KP-1 ) = TEMP 1735 END IF 1736 END IF 1737* 1738 K = K + KSTEP 1739 KC = KCNEXT 1740 GO TO 30 1741 50 CONTINUE 1742* 1743 ELSE 1744* 1745* Compute inv(A) from the factorization A = L*D*L**T. 1746* 1747* K is the main loop index, increasing from 1 to N in steps of 1748* 1 or 2, depending on the size of the diagonal blocks. 1749* 1750 NPP = N*( N+1 ) / 2 1751 K = N 1752 KC = NPP 1753 60 CONTINUE 1754* 1755* If K < 1, exit from loop. 1756* 1757 IF( K.LT.1 ) 1758 $ GO TO 80 1759* 1760 KCNEXT = KC - ( N-K+2 ) 1761 IF( IPIV( K ).GT.0 ) THEN 1762* 1763* 1 x 1 diagonal block 1764* 1765* Invert the diagonal block. 1766* 1767 AP( KC ) = ONE / AP( KC ) 1768* 1769* Compute column K of the inverse. 1770* 1771 IF( K.LT.N ) THEN 1772 CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 ) 1773 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1, 1774 $ ZERO, AP( KC+1 ), 1 ) 1775 AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 ) 1776 END IF 1777 KSTEP = 1 1778 ELSE 1779* 1780* 2 x 2 diagonal block 1781* 1782* Invert the diagonal block. 1783* 1784 T = ABS( AP( KCNEXT+1 ) ) 1785 AK = AP( KCNEXT ) / T 1786 AKP1 = AP( KC ) / T 1787 AKKP1 = AP( KCNEXT+1 ) / T 1788 D = T*( AK*AKP1-ONE ) 1789 AP( KCNEXT ) = AKP1 / D 1790 AP( KC ) = AK / D 1791 AP( KCNEXT+1 ) = -AKKP1 / D 1792* 1793* Compute columns K-1 and K of the inverse. 1794* 1795 IF( K.LT.N ) THEN 1796 CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 ) 1797 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1, 1798 $ ZERO, AP( KC+1 ), 1 ) 1799 AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 ) 1800 AP( KCNEXT+1 ) = AP( KCNEXT+1 ) - 1801 $ DDOT( N-K, AP( KC+1 ), 1, 1802 $ AP( KCNEXT+2 ), 1 ) 1803 CALL DCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 ) 1804 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1, 1805 $ ZERO, AP( KCNEXT+2 ), 1 ) 1806 AP( KCNEXT ) = AP( KCNEXT ) - 1807 $ DDOT( N-K, WORK, 1, AP( KCNEXT+2 ), 1 ) 1808 END IF 1809 KSTEP = 2 1810 KCNEXT = KCNEXT - ( N-K+3 ) 1811 END IF 1812* 1813 KP = ABS( IPIV( K ) ) 1814 IF( KP.NE.K ) THEN 1815* 1816* Interchange rows and columns K and KP in the trailing 1817* submatrix A(k-1:n,k-1:n) 1818* 1819 KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1 1820 IF( KP.LT.N ) 1821 $ CALL DSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 ) 1822 KX = KC + KP - K 1823 DO 70 J = K + 1, KP - 1 1824 KX = KX + N - J + 1 1825 TEMP = AP( KC+J-K ) 1826 AP( KC+J-K ) = AP( KX ) 1827 AP( KX ) = TEMP 1828 70 CONTINUE 1829 TEMP = AP( KC ) 1830 AP( KC ) = AP( KPC ) 1831 AP( KPC ) = TEMP 1832 IF( KSTEP.EQ.2 ) THEN 1833 TEMP = AP( KC-N+K-1 ) 1834 AP( KC-N+K-1 ) = AP( KC-N+KP-1 ) 1835 AP( KC-N+KP-1 ) = TEMP 1836 END IF 1837 END IF 1838* 1839 K = K - KSTEP 1840 KC = KCNEXT 1841 GO TO 60 1842 80 CONTINUE 1843 END IF 1844* 1845 RETURN 1846* 1847* End of DSPTRI 1848* 1849 END 1850 SUBROUTINE DSPCON( UPLO, N, AP, IPIV, ANORM, RCOND, WORK, IWORK, 1851 $ INFO ) 1852* 1853* -- LAPACK computational routine (version 3.4.0) -- 1854* -- LAPACK is a software package provided by Univ. of Tennessee, -- 1855* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 1856* November 2011 1857* 1858* .. Scalar Arguments .. 1859 CHARACTER UPLO 1860 INTEGER INFO, N 1861 DOUBLE PRECISION ANORM, RCOND 1862* .. 1863* .. Array Arguments .. 1864 INTEGER IPIV( * ), IWORK( * ) 1865 DOUBLE PRECISION AP( * ), WORK( * ) 1866* .. 1867* 1868* ===================================================================== 1869* 1870* .. Parameters .. 1871 DOUBLE PRECISION ONE, ZERO 1872 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 1873* .. 1874* .. Local Scalars .. 1875 LOGICAL UPPER 1876 INTEGER I, IP, KASE 1877 DOUBLE PRECISION AINVNM 1878* .. 1879* .. Local Arrays .. 1880 INTEGER ISAVE( 3 ) 1881* .. 1882* .. External Functions .. 1883 LOGICAL LSAME 1884 EXTERNAL LSAME 1885* .. 1886* .. External Subroutines .. 1887 EXTERNAL DLACN2, DSPTRS, XERBLA 1888* .. 1889* .. Executable Statements .. 1890* 1891* Test the input parameters. 1892* 1893 INFO = 0 1894 UPPER = LSAME( UPLO, 'U' ) 1895 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 1896 INFO = -1 1897 ELSE IF( N.LT.0 ) THEN 1898 INFO = -2 1899 ELSE IF( ANORM.LT.ZERO ) THEN 1900 INFO = -5 1901 END IF 1902 IF( INFO.NE.0 ) THEN 1903 CALL XERBLA( 'DSPCON', -INFO ) 1904 RETURN 1905 END IF 1906* 1907* Quick return if possible 1908* 1909 RCOND = ZERO 1910 IF( N.EQ.0 ) THEN 1911 RCOND = ONE 1912 RETURN 1913 ELSE IF( ANORM.LE.ZERO ) THEN 1914 RETURN 1915 END IF 1916* 1917* Check that the diagonal matrix D is nonsingular. 1918* 1919 IF( UPPER ) THEN 1920* 1921* Upper triangular storage: examine D from bottom to top 1922* 1923 IP = N*( N+1 ) / 2 1924 DO 10 I = N, 1, -1 1925 IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO ) 1926 $ RETURN 1927 IP = IP - I 1928 10 CONTINUE 1929 ELSE 1930* 1931* Lower triangular storage: examine D from top to bottom. 1932* 1933 IP = 1 1934 DO 20 I = 1, N 1935 IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO ) 1936 $ RETURN 1937 IP = IP + N - I + 1 1938 20 CONTINUE 1939 END IF 1940* 1941* Estimate the 1-norm of the inverse. 1942* 1943 KASE = 0 1944 30 CONTINUE 1945 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 1946 IF( KASE.NE.0 ) THEN 1947* 1948* Multiply by inv(L*D*L**T) or inv(U*D*U**T). 1949* 1950 CALL DSPTRS( UPLO, N, 1, AP, IPIV, WORK, N, INFO ) 1951 GO TO 30 1952 END IF 1953* 1954* Compute the estimate of the reciprocal condition number. 1955* 1956 IF( AINVNM.NE.ZERO ) 1957 $ RCOND = ( ONE / AINVNM ) / ANORM 1958* 1959 RETURN 1960* 1961* End of DSPCON 1962* 1963 END 1964#endif 1965 real*8 function ddot_128(n,da,inca,db,incb) 1966! 1967! Purpose: accumulate in quadruple precision 1968! to avoid round-off errors for big vectors 1969! hjaaj July 2012 1970! 1971 implicit none 1972 real*8 da(*), db(*), ddot 1973 integer n, inca, incb, i 1974 ! only real*16 for debugging (e.g. PGI compilers do not allow real*16) 1975 ! real*16 ddot_value 1976 real*8 ddot_value 1977 1978 if (inca .ne. 1 .or. incb .ne. 1) then 1979 ddot_128 = ddot(n,da,inca,db,incb) 1980 return 1981 end if 1982 ddot_value = 0 1983 do i = 1,n 1984 ddot_value = ddot_value + da(i)*db(i) 1985 end do 1986 ddot_128 = ddot_value 1987 return 1988 end 1989C --- end of linextra.F --- 1990