1c program qqq 2c with diis implementation 3 implicit none 4 integer nu, nv, ngr ,icr,icl,i,dd 5 real*8 r1, r2, tau, solperm,t, tol,lambd,dr 6c parameters(nv number of solute sites nu number of solute sites ngr number of grid points) 7c solperm solvent relative permittivity 8c tau(aa-1) separation between short and long range functions in 1/angstroms 9c icr type of combination rule: 1-arithmetic; 2-geometric 10c t(k) temperature in kelvins 11c tol residue tolerance 12c lamd mixing parameter 13c icl closure type: 1-hnc; 2-kh 14c tau debye parameter for spliting potnetial into short and long part 15c solperm dielectric permitivity of the solvent 16c dd number of diis basis vectors 17c ngr=4096 18c nu=12 19c nv=4 20c tau=1 21c solperm=1 22c icr=1 23c t=298 24c tol=1e-5 25c lambd=0.5 26c icl=2 27 open(3, file='uvpar.data', status='old') 28 read(3,102) nu,nv,icr,icl,ngr,tau,solperm,t,tol,lambd,dd 29 print*, nu,nv,icr,icl,ngr,tau,solperm,t,tol,lambd,dd 30 close(3) 31c parameters(r1 first grid point r2 second grid point) 32c reading of r1 and r2 33 open(3, file='full.data', status='old') 34 read(3,104) r1,r2 35 close(3) 36104 format(1x,f12.7) 37 dr=r2-r1 38! print*, r1,r2,dr 39 call risminput(nv,nu,ngr,dr,r1,tau, 40 * solperm,icr,t,tol,lambd,icl,dd) 41 102 format(4(1x,i4),1x,i8,5(2x,e9.4),1x,i4) 42 stop 43 end 44cccc 45 subroutine risminput(nv,nu,ngr,dr,r1,tau, 46 * solperm,icr,t,tol,lambd,icl,dd) 47 implicit none 48 integer nu, isu(1:nu), mu(1:nu) 49 real*8 wu(1:nu,1:nu,1:ngr),sgvv(1:nv) 50 real*8 xu(1:nu),yu(1:nu),zu(1:nu),sigu(1:nu) 51 real*8 epsiu(1:nu),qqu(1:nu) 52 integer nv, ngr, isv(1:nv), mv(1:nv) 53 real*8 den(1:nv),xv(1:nv),yv(1:nv),zv(1:nv),sigv(1:nv) 54 real*8 epsiv(1:nv),qqv(1:nv) 55 real*8 sigvv(1:nv), epsvv(1:nv),qvv(1:nv) 56 integer i, icr, icl,nvv,ims(1:nv),nd,k,dd 57 real*8 r1,dr, dk, rgrid(1:ngr),kgrid(1:ngr),pi 58 character (4) tv(1:nv), tu(1:nu) 59 real*8 tau,solperm,t, tol,lambd,kb 60 integer j1,j 61 real*8 c3(1:nu,1:nv,1:ngr),hu(1:nu,1:nv,1:ngr) 62 real*8 gam(1:nu,1:nv,1:ngr) 63c 64c solvent input 65c tv type of cite isv type of species mv multiplicity den density (1/a3) xvyvzv[xyz](a) 66c sigv sigma(a) epsiv epsilon(kj/mol) qqv charge(e) 67c number of lines in solvt..data should be equal nv 68 open(3, file='solvent3.data', status='old') 69 do i=1,nv 70 read(3,103)tv(i),isv(i),mv(i),den(i), 71 * xv(i),yv(i),zv(i),sigv(i), epsiv(i),qqv(i) 72! print *, tv(i),isv(i),mv(i),den(i), 73! * xv(i),yv(i),zv(i),sigv(i), epsiv(i),qqv(i) 74 enddo 75 close(3) 76 103 format(a4,2x,i4,2x,i4,7(2x,e10.4)) 77c 78c sorting 79 call sort(nv,tv,isv,mv,ims,nvv) 80! print*, (ims(i), i=1,nv) 81c calculations interaction parameters for reduced matrix 82 call vpot(nv,ims,nvv,sigv,epsiv,qqv,sgvv,epsvv,qvv) 83! do i=1,nvv 84! print*, sgvv(i),epsvv(i),qvv(i),i,ims(i) 85! enddo 86c 87c solute input 88c tu type of cite isu type of species mu multiplicity xuyuzu[xyz](a) 89c sigu sigma(a) epsiu epsilon(kj/mol) qqu charge(e) 90 open(3, file='solute2.data', status='old') 91 do i=1,nu 92 read(3,102)tu(i),isu(i),mu(i), 93 * xu(i),yu(i),zu(i),sigu(i),epsiu(i),qqu(i) 94 print *,tu(i),isu(i),mu(i), 95 * xu(i),yu(i),zu(i),sigu(i),epsiu(i),qqu(i) 96 enddo 97 close(3) 98 102 format(a4,2x,i4,2x,i4,6(2x,e10.4)) 99c 100c creat rgrid rgrd(i) creat kgrid kgrid(i) 101 pi=2*asin(1.0) 102c the zero point is excluded 103 r1=dr 104 dk=pi/(ngr+1)/dr 105 do i=1,ngr 106 rgrid(i)=r1+dr*(i-1) 107 kgrid(i)=i*dk 108 enddo 109c 110c wu(i,j,k) intramolecular solute function 111 call wucreat(nu,ngr,kgrid,tu,isu,xu,yu,zu,wu) 112! do k=1,ngr,40 113! print*, (wu(1,i,k), i=1,9) 114! enddo 115c 116 call rism(nu,nv,nvv,ngr,icl,kgrid,rgrid,tv,isv,den,mv,xv,yv,zv, 117 * icr,ims,tau,solperm,sgvv,epsvv,qvv,sigu,epsiu,qqu,t,lambd, 118 * tol,wu,dd) 119 return 120 end subroutine 121cccc 122c 123c subroutine for rearrangement of the solvent data which excludes the same sites 124c defined by paramters tv and mv 125c 126 subroutine sort(nv,tv,isv,mv,ims,nvv) 127 implicit none 128 integer nv, isv(1:nv), mv(1:nv),ims(1:nv),mvt(1:nv) 129 integer i,i1,j1,j2,icc,imm(1:nv),max,nvv 130 character (4) tv(1:nv) 131c 132c initialization of connectivity vector 133 do i=1,nv 134 ims(i)=i 135 imm(i)=1 136 mvt(i)=mv(i) 137 enddo 138c sorting 139 do j1=1,nv-1 140 do j2=j1+1,nv 141 if((tv(j1).eq.tv(j2)).and.(isv(j1).eq.isv(j2)) 142 * .and.(mv(j2).ne.1)) then 143c replace current value of ims by the first found 144 ims(j2)=ims(j1) 145c change indicator of replacement 146 imm(j2)=2 147c change multiplicity 148 mv(j2)=1 149c shift all others except replaced by 1 150 if(j2.ne.nv) then 151 do i1=j2+1,nv 152 if(imm(i1).eq.1) then 153 ims(i1)=ims(i1)-1 154 else 155 ims(i1)=ims(i1) 156 endif 157 enddo 158 endif 159 else 160 ims(j2)=ims(j2) 161 endif 162! print*, j2,ims(j2),mv(j2),imm(j2) 163 enddo 164 enddo 165c evaluation of maximal value of different sites 166 max=ims(1) 167 do i=2,nv 168 if(ims(i).ge.max) then 169 max=ims(i) 170 else 171 max=max 172 endif 173 enddo 174 nvv=max 175c replace changed mv by the initial one 176 do i=1,nv 177 mv(i)=mvt(i) 178 enddo 179 return 180 end subroutine 181c 182c subroutine for calculation potential parameters for reduced matrix 183c 184 subroutine vpot(nv,ims,nvv,sigv,epsiv,qqv,sgvv,epsvv,qvv) 185 implicit none 186 integer nv,nvv 187 real*8 sgvv(1:nvv),epsvv(1:nvv),qvv(1:nvv) 188 real*8 sigv(1:nv),epsiv(1:nv),qqv(1:nv) 189 integer i,ims(1:nv),j1,j2,icr 190 do i=1,nv 191 sgvv(ims(i))=sigv(i) 192 epsvv(ims(i))=epsiv(i) 193 qvv(ims(i))=qqv(i) 194! print*, sgvv(ims(i)),sigv(i),i,ims(i) 195 enddo 196 return 197 end subroutine 198c 199c prepration of wu(i,j,k) intramoplecular solute matrix 200c 201 subroutine wucreat(nu,ngr,kgrid,tu,isu,xu,yu,zu,wu) 202 implicit none 203 integer nu, ngr, isu(1:nu), mu(1:nu) 204 real*8 wu(1:nu,1:nu,1:ngr), dist(1:nu,1:nu) 205 real*8 xu(1:nu),yu(1:nu),zu(1:nu) 206 integer i, j1,j2,icr 207 real*8 kgrid(1:ngr),sik 208 character (4) tu(1:nu) 209c dis(i,j) intramolecular distances 210 do i=1,ngr 211 do j1=1,nu 212 do j2=1,nu 213 dist(j1,j2)=(((xu(j1)-xu(j2))**2+(yu(j1)-yu(j2))**2+ 214 * (zu(j1)-zu(j2))**2))**(1.0/2) 215 wu(j1,j2,i)=sik(kgrid(i)*dist(j1,j2)) 216 enddo 217 enddo 218 enddo 219 return 220 end subroutine 221c 222 real*8 function sik(x) 223 real*8 x 224 if(x.lt.1d-22)then 225 sik=1 226 else 227 sik=sin(x)/x 228 endif 229 return 230 end 231c 232c creat susceptibility of solvent 233c 234 subroutine chicreat(nd,nv,ngr,ims,nvv,kgrid,tv, 235 * isv,den,mv,xv,yv,zv,chi) 236 implicit none 237 real*8 wv(1:nvv,1:nvv,1:ngr), chi(1:nvv,1:nvv,1:ngr),hd(1:ngr) 238 integer nv, ngr, isv(1:nv), mv(1:nv) 239 real*8 den(1:nv),xv(1:nv),yv(1:nv),zv(1:nv),sigv(1:nv) 240 integer i,j,j1,j2,icr,nvv,ims(1:nv),nd 241 real*8 r1, r2, dr, dk, rgrid(1:ngr),kgrid(1:ngr),pi 242 real*8 hsol(1:nd,1:ngr), rr(1:ngr) 243 real*8 h1(1:nvv,1:nvv,1:ngr),h2(1:ngr),h3(1:nvv,1:nvv,1:ngr) 244 character (4) tv(1:nv) 245 pi=2*asin(1.0) 246 247c 248c reading of solvent rdfs 249 open(3, file='full.data', status='old') 250 do i=1,ngr 251 read(3,104) rr(i), (hsol(j,i), j=1,nd) 252 enddo 253 close(3) 254104 format(11(1x,f12.7)) 255 dr=rr(2)-rr(1) 256! do i=1,ngr,40 257! print*, (hsol(j,i), j=1,9) 258! enddo 259c 260c sin fft transform of rdfs with proper arrangement in array 261 do j=1,nvv 262 do j1=j,nvv 263 do i=2,ngr 264 h1(j,j1,i)=(hsol(nvv*(j-1)-j*(j-1)/2+j1,i)-1)*rr(i) 265 h2(i)=h1(j,j1,i) 266 enddo 267 call sinft(h2,ngr) 268c normalization of sin-fft with excluding the zeropoint (x=0) 269 do i=1,ngr-1 270 h3(j,j1,i)=h2(i+1)/kgrid(i) 271 enddo 272 h3(j,j1,ngr)=h3(j,j1,ngr-1) 273 enddo 274 enddo 275c symmetric indeces 276 do j=1,nvv 277 do j1=j,nvv 278 do i=1,ngr 279 h3(j1,j,i)=h3(j,j1,i) 280 enddo 281 enddo 282 enddo 283! do i=1,ngr,40 284! print*, (h3(j,3,i), j=1,4) 285! enddo 286c 287c wv(i,j,k) intramolecular solvent function 288c 289 call wcreat(nv,ngr,nvv,ims,kgrid,tv,isv,xv,yv,zv,wv) 290! do i=1,ngr,40 291! print*, (wv(1,j1,i), j1=1,nvv) 292! enddo 293c 294c chi(i,j,k) solvent susceptibility 295c 296 dr=rr(2)-rr(1) 297 do i=1,ngr 298 do j1=1,nv 299 do j2=1,nv 300 chi(ims(j1),ims(j2),i)=wv(ims(j1),ims(j2),i) 301 * +4*pi*dr*mv(j1)*mv(j2)*den(j1)*h3(ims(j1),ims(j2),i) 302 enddo 303 enddo 304 enddo 305! do i=1,ngr,40 306! print*, (chi(2,j,i), j=1,4) 307! enddo 308 return 309 end subroutine 310c 311c caculation of reduced solvent matrix 312c 313 subroutine wcreat(nv,ngr,nvv,ims,kgrid,tv,isv,xv,yv,zv,wv) 314 implicit none 315 real*8 wv(1:nvv,1:nvv,1:ngr),dist(1:nv,1:nv) 316 real*8 ws(1:nv,1:nvv,1:ngr) 317 integer nv, ngr, isv(1:nv),mv(1:nv) 318 real*8 den(1:nv),xv(1:nv),yv(1:nv),zv(1:nv),sigv(1:nv) 319 integer i, j1,j2,nvv,ims(1:nv),t 320 real*8 kgrid(1:ngr),pi,sik,co(1:nv,1:nv,1:ngr) 321 character (4) tv(1:nv) 322 pi=2*asin(1.0) 323c 324c initialization of wv(i,j,k) & calculating of the 11 element 325c 326 do i=1, ngr 327 do j1=1,nv 328 do j2=1,nv 329 wv(ims(j1),ims(j2),i)=0 330 ws(j1,ims(j2),i)=0 331 enddo 332 enddo 333 enddo 334c dist(i,j) distance between i and j sites 335c co(i,j,k) sik between i and j sites 336c reducing number of coloumns 337 do i=1,ngr 338 do j1=1,nv 339 do j2=1,nv 340c calculations all other columns 341c if they belong to the same molecule 342 if(isv(j1).eq.isv(j2)) then 343 dist(j1,j2)=(((xv(j1)-xv(j2))**2+ 344 * (yv(j1)-yv(j2))**2+(zv(j1)-zv(j2))**2))**(1.0/2) 345 co(j1,j2,i)=sik(kgrid(i)*dist(j1,j2)) 346 ws(j1,ims(j2),i)= ws(j1,ims(j2),i)+co(j1,j2,i) 347 endif 348 enddo 349 enddo 350 enddo 351c reducing the number of rows 352 do i=1,ngr 353 do j2=1,nvv 354 do j1=1,nv 355 wv(ims(j1),j2,i)=wv(ims(j1),j2,i)+ws(j1,j2,i) 356 enddo 357 enddo 358 enddo 359! do i=1,ngr,40 360! print*, (wv(1,j1,i), j1=1,nvv) 361! enddo 362 return 363 end subroutine 364c 365c creat solute-solvent potentials 366c 367 subroutine potcreat(nvv,nu,ngr,icr,kgrid,rgrid,tau,solperm, 368 * sgvv,epsvv,qvv,sigu,epsiu,qqu,plj,ul) 369 implicit none 370 integer nu,nv, ngr,nvv 371 real*8 sigu(1:nu),epsiu(1:nu),qqu(1:nu) 372 real*8 sgvv(1:nvv),epsvv(1:nvv),qvv(1:nvv) 373 integer i, j1,j2,icr 374 real*8 r1, r2, dr, dk, rgrid(1:ngr),kgrid(1:ngr),pi 375 real*8 sigvv(1:nu,1:nvv), epsivv(1:nu,1:nvv) 376 real*8 ul(1:nu,1:nvv,1:ngr),plj(1:nu,1:nvv,1:ngr) 377 real*8 nav,tau, echar,solperm,uthre,ck 378c constants for potential 379 pi=2*asin(1.0) 380 nav=6.02214179e+23 381 echar=4.803e-10 382 ck=nav*echar**2/solperm*0.01 383 uthre=1000 384c 385c calculation lj parameters 386 call combrule(nu,nvv,icr,sgvv,sigu,epsiu,epsvv,sigvv,epsivv) 387c 388c plj(j1,j2,r) lj potential+short part from coulomb 389c ul(j1,j2,k) is the fourier trasnform of long range interactions 390 do i=1,ngr 391 do j1=1,nu 392 do j2=1,nvv 393 ul(j1,j2,i)=4*pi*ck*qqu(j1)*qvv(j2)*exp(-kgrid(i)**2/4/tau**2) 394 ul(j1,j2,i)=ul(j1,j2,i)/kgrid(i)**2 395 plj(j1,j2,i)=4*epsivv(j1,j2)* 396 * ((sigvv(j1,j2)/rgrid(i))**(12)-(sigvv(j1,j2)/rgrid(i))**(6)) 397 plj(j1,j2,i)=plj(j1,j2,i)+ck 398 * *qqu(j1)*qvv(j2)/rgrid(i)*(1-erf(tau*rgrid(i))) 399 if (plj(j1,j2,i).ge. uthre) then 400 plj(j1,j2,i)=uthre 401 endif 402 enddo 403 enddo 404 enddo 405! do i=1,40 406! print*, (plj(1,j1,i), j1=1,4) 407! enddo 408 return 409 end subroutine 410c 411c calculations of parameters for solute-solvent interactions 412c 413 subroutine combrule(nu,nvv,icr,sgvv,sigu, 414 * epsiu,epsvv,sigvv,epsivv) 415 implicit none 416 integer nu, nv, icr,nvv 417 real*8 sigu(1:nu), epsiu(1:nu), sgvv(1:nvv), epsvv(1:nvv) 418 real*8 sigvv(1:nu,1:nvv), epsivv(1:nu,1:nvv) 419 integer i, j1,j2 420 do j1=1,nu 421 do j2=1,nvv 422 epsivv(j1,j2)=(epsiu(j1)*epsvv(j2))**(1.0/2) 423 if (icr.eq.1) then 424 sigvv(j1,j2)=(sigu(j1)+sgvv(j2))/2 425 else 426 sigvv(j1,j2)=(sigu(j1)*sgvv(j2))**(1.0/2) 427 endif 428! print *, sigvv(j1,j2), epsivv(j1,j2), j1,j2 429 enddo 430 enddo 431 return 432 end subroutine 433c 434c site-site ornstein-zernike in k-space 435c 436 subroutine ssoz(nvv,nu,ngr,nv,ims,mv,wu,c3,chi,hu) 437 implicit none 438 integer nu, nvv, nv, ngr, i,j1,j,k1,ims(1:nv),mv(1:nv) 439 real*8 rgrid(1:ngr),kgrid(1:ngr) 440 real*8 pi,ct 441 real*8 c3(1:nu,1:nvv,1:ngr),hu(1:nu,1:nvv,1:ngr) 442 real*8 ctt(1:nu,1:nvv,1:ngr),h1(1:nu,1:nvv,1:ngr) 443 real*8 wu(1:nu,1:nu,1:ngr), chi(1:nvv,1:nvv,1:ngr) 444c calculations of right product of matrix via summation over solvent sites 445 do i=1,ngr 446 do j1=1,nu 447 do j=1,nvv 448 ct=0 449 do k1=1,nvv 450 ct=c3(j1,k1,i)*chi(k1,j,i)+ct 451 enddo 452 ctt(j1,j,i)=ct 453 enddo 454 enddo 455 enddo 456! do i=1,ngr,40 457! print*, (ctt(1,j,i), j=1,4) 458! enddo 459c calculations of left product of matrix via summation over solute sites 460 do i=1,ngr 461 do j1=1,nu 462 do j=1,nvv 463 ct=0 464 do k1=1,nu 465 ct=wu(j1,k1,i)*ctt(k1,j,i)+ct 466 enddo 467 hu(j1,j,i)=ct 468 h1(j1,j,i)=hu(j1,j,i) 469 enddo 470 enddo 471 enddo 472 do i=1,ngr 473 do j1=1,nu 474 do j=1,nv 475 hu(j1,ims(j),i)=h1(j1,ims(j),i)/mv(j) 476 enddo 477 enddo 478 enddo 479! do i=1,ngr,40 480! print*, (hu(1,j,i),h1(1,j,i), j=1,4) 481! enddo 482 return 483 end subroutine 484c 485c subroutine for closure 486c 487 subroutine clos(nvv,nu,ngr,t,plj,cr,gam,icl) 488 implicit none 489 integer nu, nvv, ngr, i,j1,j,icl 490 real*8 plj(1:nu,1:nvv,1:ngr) 491 real*8 t, kb,pi,lexp 492 real*8 cr(1:nu,1:nvv,1:ngr),gam(1:nu,1:nvv,1:ngr) 493 kb=8.13441e-3 494 do i=1,ngr 495 do j=1,nu 496 do j1=1,nvv 497 if(icl.eq.1) then 498 cr(j,j1,i)=exp(-1.0/kb/t*plj(j,j1,i)+gam(j,j1,i)) 499 * -gam(j,j1,i)-1 500 endif 501 if(icl.eq.2) then 502 cr(j,j1,i)=lexp(-1.0/kb/t*plj(j,j1,i)+gam(j,j1,i)) 503 * -gam(j,j1,i)-1 504 endif 505 enddo 506 enddo 507 enddo 508 return 509 end subroutine 510c 511 real*8 function lexp(x) 512 real*8 x 513 if(x.le.0.0)then 514 lexp=exp(x) 515 else 516 lexp=1+x 517 endif 518 return 519 end 520c 521ccc 522c 523 subroutine rism(nu,nv,nvv,ngr,icl,kgrid,rgrid,tv,isv, 524 * den,mv,xv,yv,zv, 525 * icr,ims,tau,solperm,sgvv,epsvv,qvv,sigu,epsiu,qqu,t, 526 * lambd,tol,wu,dd) 527 implicit none 528c parameters 529 integer nu, nv,nd,nvv, ngr, i,j1,j,k1,icl,k 530 real*8 rgrid(1:ngr),kgrid(1:ngr) 531 real*8 pi,dk,t,tau,solperm,tol,lambd,del0,kb,dr 532c solute 533 real*8 wu(1:nu,1:nu,1:ngr), chi(1:nvv,1:nvv,1:ngr) 534 real*8 sigu(1:nu),epsiu(1:nu),qqu(1:nu) 535c solvent 536 integer icr,isv(1:nv),mv(1:nv),ims(1:nv) 537 real*8 den(1:nv),xv(1:nv),yv(1:nv),zv(1:nv) 538 real*8 wv(1:nvv,1:nvv,1:ngr) 539 real*8 sgvv(1:nvv),epsvv(1:nvv),qvv(1:nvv) 540 character (4) tv(1:nv) 541c potential 542 real*8 plj(1:nu,1:nvv,1:ngr), ul(1:nu,1:nvv,1:ngr) 543c functions 544 real*8 cr(1:nu,1:nvv,1:ngr),c2(1:ngr),tt(1:nu,1:nvv) 545 real*8 ht(1:ngr),c3(1:nu,1:nvv,1:ngr),cf3(1:nu,1:nvv,1:ngr) 546 real*8 cold(1:nu,1:nvv,1:ngr),cnew(1:nu,1:nvv,1:ngr) 547 real*8 gfold(1:nu,1:nvv,1:ngr),hu(1:nu,1:nvv,1:ngr) 548 real*8 gold(1:nu,1:nvv,1:ngr), gnew(1:nu,1:nvv,1:ngr) 549 real*8 hold(1:nu,1:nvv,1:ngr), hnew(1:nu,1:nvv,1:ngr) 550 real*8 del(1:nu,1:nvv),mmaxi,mmax(1:nu),mmmax 551c diis functions and counts 552 real*8 ggo(1:dd,1:nu,1:nvv,1:ngr),dgg(1:dd,1:nu,1:nvv,1:ngr) 553 integer kd,dd 554c 555 pi=2*asin(1.0) 556c nd total number of different rdfs 557 nd=(nvv+1)*nvv/2 558c chi(i,j,k) solvent susceptibility 559c 560 call chicreat(nd,nv,ngr,ims,nvv,kgrid,tv,isv,den,mv,xv,yv,zv,chi) 561! do i=1,ngr,40 562! print*, (chi(2,j,i), j=1,4) 563! enddo 564c 565c plj(j1,j2,r) lj potential+short part from coulomb 566c ul(j1,j2,k) is the fourier trasnform of long range interactions 567c 568 call potcreat(nvv,nu,ngr,icr,kgrid,rgrid,tau,solperm, 569 * sgvv,epsvv,qvv,sigu,epsiu,qqu,plj,ul) 570! do k=1,40 571! print*, (plj(i,1,k), i=1,9) 572! enddo 573cc initial value of c_short(r) 574 kb=8.13441e-3 575 do i=1,ngr 576 do j=1,nu 577 do j1=1,nvv 578 gold(j,j1,i)=0 579 cold(j,j1,i)=0 580 enddo 581 enddo 582 enddo 583c starting counts for diis (kd) and for iterations k1 584 k1=1 585 kd=1 586 1 continue 587c calculating c_new(r) and h_new(r) from closure 588 call clos(nvv,nu,ngr,t,plj,cnew,gold,icl) 589 do j=1,nu 590 do j1=1,nvv 591 do i=1,ngr 592 hold(j,j1,i)=gold(j,j1,i)+cnew(j,j1,i) 593 enddo 594 enddo 595 enddo 596! do i=1,ngr,40 597! print*, rgrid(i),(cnew(1,j,i), j=1,4) 598! enddo 599c fast fourier transform 600 dr=rgrid(2)-rgrid(1) 601 do j=1,nu 602 do j1=1,nvv 603 c2(1)=0 604 do i=1,ngr-1 605 c2(i+1)=cnew(j,j1,i)*rgrid(i) 606 enddo 607 call sinft(c2,ngr) 608c normalization of sin-fft with excluding the zeropoint (x=0) 609 do i=1,ngr-1 610 cf3(j,j1,i)=4*pi*c2(i+1)/kgrid(i)*dr 611 enddo 612 cf3(j,j1,ngr)=cf3(j,j1,ngr-1) 613 enddo 614 enddo 615! do i=1,40 616! print*, kgrid(i), (cf3(1,j1,i), j1=1,4) 617! enddo 618 dk=kgrid(2)-kgrid(1) 619 pi=2*asin(1.0) 620c adding the long-range potential to c_short 621c 622 do i=1,ngr 623 do j=1,nu 624 do j1=1,nvv 625 cf3(j,j1,i)=cf3(j,j1,i)-1.0/kb/t*ul(j,j1,i) 626 enddo 627 enddo 628 enddo 629c calculations of fourier transform of h(i,j,k) by 630c site-site ornstein-zerinike equations 631c 632 call ssoz(nvv,nu,ngr,nv,ims,mv,wu,cf3,chi,hu) 633! do i=1,ngr,40 634! print*, (hu(1,j,i), j=1,4) 635! enddo 636c inverse fourier transform of h(i,j,k) & calculations of gamma(i,i,r) 637 do j=1,nu 638 do j1=1,nvv 639 ht(1)=0 640 do i=1,ngr-1 641 ht(i+1)=hu(j,j1,i)*kgrid(i) 642 enddo 643 call sinft(ht,ngr) 644 do i=1,ngr-1 645 hnew(j,j1,i)=ht(i+1)*dk/2/rgrid(i)/pi**2 646 gnew(j,j1,i)=hnew(j,j1,i)-cnew(j,j1,i) 647 enddo 648 gnew(j,j1,ngr)=gnew(j,j1,ngr-1) 649 hnew(j,j1,ngr)=hnew(j,j1,ngr-1) 650 enddo 651 enddo 652! do i=1,40 653! print*, (gnew(1,j,i), j=1,4) 654! enddo 655c intialization del 656 do j=1,nu 657 do j1=1,nvv 658 del(j,j1)=0 659 enddo 660 mmax(j)=0 661 enddo 662 mmaxi=0 663c evaluation of accuracy 664 do j=1,nu 665 do j1=1,nvv 666 do i=1,ngr 667 del(j,j1)=del(j,j1)+(gnew(j,j1,i)-gold(j,j1,i))**2 668 enddo 669 mmax(j)=del(j,j1)+mmax(j) 670 enddo 671 mmaxi=mmaxi+mmax(j) 672 enddo 673! do j=1,nu 674! print*, (del(j,j1), j1=1,4),mmmax 675! enddo 676c normalization 677 del0=(mmaxi/ngr/nu/nvv)**(1.0/2) 678c diis procedure 679 call diis(nu,nvv,ngr,lambd,kd,dd,gold,gnew,ggo,dgg) 680c old functions as new functions for new cycle 681 do j=1,nu 682 do j1=1,nvv 683 do i=1,ngr 684 gold(j,j1,i)=gnew(j,j1,i) 685 cold(j,j1,i)=cnew(j,j1,i) 686 enddo 687 enddo 688 enddo 689 107 format(9(1x,e10.4)) 690 k1=k1+1 691 print*, k1,kd,del0 692 if(k1.ge.700) goto 2 693 if(del0.ge.tol) goto 1 694 2 continue 695 open(3, file='g1.data', status='old') 696 do k=1,ngr 697 write(3,107)(gnew(1,i,k)+cnew(1,i,k), i=1,4) 698 enddo 699 close(3) 700 return 701 end subroutine 702c 703c subroutine for evaluation diis residues dgg and matrix coefficients ss 704c 705 subroutine diis(nu,nvv,ngr,lam,k1,dd,hold,hnew,ggo,dgg) 706 implicit none 707 real*8 ggo(1:dd,1:nu,1:nvv,1:ngr),dgg(1:dd,1:nu,1:nvv,1:ngr) 708 real*8 ss(1:dd+1,1:dd+1),hold(1:nu,1:nvv,1:ngr) 709 real*8 s0(1:dd+1),s1,lam,hnew(1:nu,1:nvv,1:ngr) 710 integer k1,dd,jd,jd1 711 integer j,j1,i,nu,nvv,ngr 712 integer ipiv(1:dd+1),info 713c initialization of overlap matrix for diis 714 do jd=1,dd 715 do jd1=1,dd 716 ss(jd,jd1)=0 717 enddo 718 s0(jd)=0 719 ss(jd,dd+1)=-1 720 ss(dd+1,jd)=-1 721 enddo 722 ss(dd+1,dd+1)=0 723 s0(dd+1)=-1 724c calculation of basis 725 do j=1,nu 726 do j1=1,nvv 727 do i=1,ngr 728 ggo(k1,j,j1,i)=hnew(j,j1,i) 729 dgg(k1,j,j1,i)=hnew(j,j1,i)-hold(j,j1,i) 730c calculation of elements of overlaping matrix 731 if(k1.eq.dd) then 732 do jd=1,dd 733 do jd1=jd,dd 734 ss(jd,jd1)=(dgg(jd,j,j1,i)*dgg(jd1,j,j1,i))/ngr+ss(jd,jd1) 735c symmerization 736 ss(jd1,jd)=ss(jd,jd1) 737 enddo 738 enddo 739 endif 740 enddo 741 enddo 742 enddo 743c solution of diis equation and change of baisis 744 if(k1.eq.dd) then 745 do jd=1,dd 746c print*, (ss(jd,jd1), jd1=1,dd) 747 enddo 748 call dgesv(dd+1,1,ss,dd+1,ipiv,s0,dd+1,info) 749 do i=1,ngr 750 do j=1,nu 751 do j1=1,nvv 752 s1=0 753c calculation sum for new solution 754 do jd=1,dd 755 s1=s1+s0(jd)*ggo(jd,j,j1,i)+lam*s0(jd)*dgg(jd,j,j1,i) 756 enddo 757c change of diis basis 758 do jd=1,dd-1 759 ggo(jd,j,j1,i)=ggo(jd+1,j,j1,i) 760 dgg(jd,j,j1,i)=dgg(jd+1,j,j1,i) 761 enddo 762 hnew(j,j1,i)=s1 763 enddo 764 enddo 765 enddo 766 endif 767c matrix singular or has illegal value see dgeesv 768 if((k1.eq.dd).and.(info.ne.0)) then 769 k1=1 770 print*, info 771 go to 1 772 endif 773c change count for diis 774 if(k1.lt.dd) then 775 k1=k1+1 776 else 777 k1=dd 778 endif 779 1 continue 780 return 781 end subroutine 782c 783c uses realft 784c calculates the sine transform of a set of n real-valued data points stored in array y(1:n). 785c the number n must be a power of 2. on exit y is replaced by its transform. this program, 786c without changes, also calculates the inverse sine transform, but in this case the output array 787c should be multiplied by 2/n. 788c 789 subroutine sinft(y,n) 790 integer n 791 real*8 y(n) 792 integer j 793 real*8 sum,y1,y2 794 double precision theta,wi,wpi,wpr,wr,wtemp 795 theta=3.141592653589793d0/dble(n) 796 wr=1.0d0 797 wi=0.0d0 798 wpr=-2.0d0*sin(0.5d0*theta)**2 799 wpi=sin(theta) 800 y(1)=0.0 801 do 11 j=1,n/2 802 wtemp=wr 803 wr=wr*wpr-wi*wpi+wr 804 wi=wi*wpr+wtemp*wpi+wi 805 y1=wi*(y(j+1)+y(n-j+1)) 806 y2=0.5*(y(j+1)-y(n-j+1)) 807 y(j+1)=y1+y2 808 y(n-j+1)=y1-y2 80911 continue 810 call realft(y,n,+1) 811 sum=0.0 812 y(1)=0.5*y(1) 813 y(2)=0.0 814 do 12 j=1,n-1,2 815 sum=sum+y(j) 816 y(j)=y(j+1) 817 y(j+1)=sum 81812 continue 819 return 820 end subroutine 821c 822cc uses four1 823c calculates the fourier transform of a set of n real-valued data points. replaces this data 824c (which is stored in array data(1:n)) by the positive frequency half of its complex fourier 825c transform. the real-valued rst and last components of the complex transform are returned 826c as elements data(1) and data(2), respectively. n must be a power of 2. this routine 827c also calculates the inverse transform of a complex data array if it is the transform of real 828c data. (result in this case must be multiplied by 2/n.) 829c 830 subroutine realft(data,n,isign) 831 integer isign,n 832 real*8 data(n) 833 integer i,i1,i2,i3,i4,n2p3 834 real*8 c1,c2,h1i,h1r,h2i,h2r,wis,wrs 835 double precision theta,wi,wpi,wpr,wr,wtemp 836 theta=3.141592653589793d0/dble(n/2) 837 c1=0.5 838 if (isign.eq.1) then 839 c2=-0.5 840 call four1(data,n/2,+1) 841 else 842 c2=0.5 843 theta=-theta 844 endif 845 wpr=-2.0d0*sin(0.5d0*theta)**2 846 wpi=sin(theta) 847 wr=1.0d0+wpr 848 wi=wpi 849 n2p3=n+3 850 do 11 i=2,n/4 851 i1=2*i-1 852 i2=i1+1 853 i3=n2p3-i2 854 i4=i3+1 855 wrs=sngl(wr) 856 wis=sngl(wi) 857 h1r=c1*(data(i1)+data(i3)) 858 h1i=c1*(data(i2)-data(i4)) 859 h2r=-c2*(data(i2)+data(i4)) 860 h2i=c2*(data(i1)-data(i3)) 861 data(i1)=h1r+wrs*h2r-wis*h2i 862 data(i2)=h1i+wrs*h2i+wis*h2r 863 data(i3)=h1r-wrs*h2r+wis*h2i 864 data(i4)=-h1i+wrs*h2i+wis*h2r 865 wtemp=wr 866 wr=wr*wpr-wi*wpi+wr 867 wi=wi*wpr+wtemp*wpi+wi 86811 continue 869 if (isign.eq.1) then 870 h1r=data(1) 871 data(1)=h1r+data(2) 872 data(2)=h1r-data(2) 873 else 874 h1r=data(1) 875 data(1)=c1*(h1r+data(2)) 876 data(2)=c1*(h1r-data(2)) 877 call four1(data,n/2,-1) 878 endif 879 return 880 end subroutine 881cc replaces data(1:2*nn) by its discrete fourier transform, if isign is input as 1; or replaces 882c data(1:2*nn) by nn times its inverse discrete fourier transform, if isign is input as -1. 883c data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn 884c must be an integer power of 2 (this is not checked for!). 885 subroutine four1(data,nn,isign) 886 integer isign,nn 887 real*8 data(2*nn) 888 integer i,istep,j,m,mmax,n 889 real*8 tempi,tempr 890 double precision theta,wi,wpi,wpr,wr,wtemp 891 n=2*nn 892 j=1 893 do 11 i=1,n,2 894 if(j.gt.i)then 895 tempr=data(j) 896 tempi=data(j+1) 897 data(j)=data(i) 898 data(j+1)=data(i+1) 899 data(i)=tempr 900 data(i+1)=tempi 901 endif 902 m=n/2 9031 if ((m.ge.2).and.(j.gt.m)) then 904 j=j-m 905 m=m/2 906 goto 1 907 endif 908 j=j+m 90911 continue 910 mmax=2 9112 if (n.gt.mmax) then 912 istep=2*mmax 913 theta=6.28318530717959d0/(isign*mmax) 914 wpr=-2.d0*sin(0.5d0*theta)**2 915 wpi=sin(theta) 916 wr=1.d0 917 wi=0.d0 918 do 13 m=1,mmax,2 919 do 12 i=m,n,istep 920 j=i+mmax 921 tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1) 922 tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j) 923 data(j)=data(i)-tempr 924 data(j+1)=data(i+1)-tempi 925 data(i)=data(i)+tempr 926 data(i+1)=data(i+1)+tempi 92712 continue 928 wtemp=wr 929 wr=wr*wpr-wi*wpi+wr 930 wi=wi*wpr+wtemp*wpi+wi 93113 continue 932 mmax=istep 933 goto 2 934 endif 935 return 936 end subroutine 937ccc lapack subroutines 938ccc 939 subroutine dgesv( n, nrhs, a, lda, ipiv, b, ldb, info ) 940* 941* -- lapack driver routine (version 3.2) -- 942* -- lapack is a software package provided by univ. of tennessee, -- 943* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 944* november 2006 945* 946* .. scalar arguments .. 947 integer info, lda, ldb, n, nrhs 948* .. 949* .. array arguments .. 950 integer ipiv( * ) 951 double precision a( lda, * ), b( ldb, * ) 952* .. 953* 954* purpose 955* ======= 956* 957* dgesv computes the solution to a real system of linear equations 958* a * x = b, 959* where a is an n-by-n matrix and x and b are n-by-nrhs matrices. 960* 961* the lu decomposition with partial pivoting and row interchanges is 962* used to factor a as 963* a = p * l * u, 964* where p is a permutation matrix, l is unit lower triangular, and u is 965* upper triangular. the factored form of a is then used to solve the 966* system of equations a * x = b. 967* 968* arguments 969* ========= 970* 971* n (input) integer 972* the number of linear equations, i.e., the order of the 973* matrix a. n >= 0. 974* 975* nrhs (input) integer 976* the number of right hand sides, i.e., the number of columns 977* of the matrix b. nrhs >= 0. 978* 979* a (input/output) double precision array, dimension (lda,n) 980* on entry, the n-by-n coefficient matrix a. 981* on exit, the factors l and u from the factorization 982* a = p*l*u; the unit diagonal elements of l are not stored. 983* 984* lda (input) integer 985* the leading dimension of the array a. lda >= max(1,n). 986* 987* ipiv (output) integer array, dimension (n) 988* the pivot indices that define the permutation matrix p; 989* row i of the matrix was interchanged with row ipiv(i). 990* 991* b (input/output) double precision array, dimension (ldb,nrhs) 992* on entry, the n-by-nrhs matrix of right hand side matrix b. 993* on exit, if info = 0, the n-by-nrhs solution matrix x. 994* 995* ldb (input) integer 996* the leading dimension of the array b. ldb >= max(1,n). 997* 998* info (output) integer 999* = 0: successful exit 1000* < 0: if info = -i, the i-th argument had an illegal value 1001* > 0: if info = i, u(i,i) is exactly zero. the factorization 1002* has been completed, but the factor u is exactly 1003* singular, so the solution could not be computed. 1004* 1005* ===================================================================== 1006* 1007* .. external subroutines .. 1008 external dgetrf, dgetrs, xerbla 1009* .. 1010* .. intrinsic functions .. 1011 intrinsic max 1012* .. 1013* .. executable statements .. 1014* 1015* test the input parameters. 1016* 1017 info = 0 1018 if( n.lt.0 ) then 1019 info = -1 1020 else if( nrhs.lt.0 ) then 1021 info = -2 1022 else if( lda.lt.max( 1, n ) ) then 1023 info = -4 1024 else if( ldb.lt.max( 1, n ) ) then 1025 info = -7 1026 end if 1027 if( info.ne.0 ) then 1028 call xerbla( 'dgesv ', -info ) 1029 return 1030 end if 1031* 1032* compute the lu factorization of a. 1033* 1034 call dgetrf( n, n, a, lda, ipiv, info ) 1035 if( info.eq.0 ) then 1036* 1037* solve the system a*x = b, overwriting b with x. 1038* 1039 call dgetrs( 'no transpose', n, nrhs, a, lda, ipiv, b, ldb, 1040 $ info ) 1041 end if 1042 return 1043* 1044* end of dgesv 1045* 1046 end 1047ccc 1048 subroutine dgetrf( m, n, a, lda, ipiv, info ) 1049* 1050* -- lapack routine (version 3.2) -- 1051* -- lapack is a software package provided by univ. of tennessee, -- 1052* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 1053* november 2006 1054* 1055* .. scalar arguments .. 1056 integer info, lda, m, n 1057* .. 1058* .. array arguments .. 1059 integer ipiv( * ) 1060 double precision a( lda, * ) 1061* .. 1062* 1063* purpose 1064* ======= 1065* 1066* dgetrf computes an lu factorization of a general m-by-n matrix a 1067* using partial pivoting with row interchanges. 1068* 1069* the factorization has the form 1070* a = p * l * u 1071* where p is a permutation matrix, l is lower triangular with unit 1072* diagonal elements (lower trapezoidal if m > n), and u is upper 1073* triangular (upper trapezoidal if m < n). 1074* 1075* this is the right-looking level 3 blas version of the algorithm. 1076* 1077* arguments 1078* ========= 1079* 1080* m (input) integer 1081* the number of rows of the matrix a. m >= 0. 1082* 1083* n (input) integer 1084* the number of columns of the matrix a. n >= 0. 1085* 1086* a (input/output) double precision array, dimension (lda,n) 1087* on entry, the m-by-n matrix to be factored. 1088* on exit, the factors l and u from the factorization 1089* a = p*l*u; the unit diagonal elements of l are not stored. 1090* 1091* lda (input) integer 1092* the leading dimension of the array a. lda >= max(1,m). 1093* 1094* ipiv (output) integer array, dimension (min(m,n)) 1095* the pivot indices; for 1 <= i <= min(m,n), row i of the 1096* matrix was interchanged with row ipiv(i). 1097* 1098* info (output) integer 1099* = 0: successful exit 1100* < 0: if info = -i, the i-th argument had an illegal value 1101* > 0: if info = i, u(i,i) is exactly zero. the factorization 1102* has been completed, but the factor u is exactly 1103* singular, and division by zero will occur if it is used 1104* to solve a system of equations. 1105* 1106* ===================================================================== 1107* 1108* .. parameters .. 1109 double precision one 1110 parameter ( one = 1.0d+0 ) 1111* .. 1112* .. local scalars .. 1113 integer i, iinfo, j, jb, nb 1114* .. 1115* .. external subroutines .. 1116 external dgemm, dgetf2, dlaswp, dtrsm, xerbla 1117* .. 1118* .. external functions .. 1119 integer ilaenv 1120 external ilaenv 1121* .. 1122* .. intrinsic functions .. 1123 intrinsic max, min 1124* .. 1125* .. executable statements .. 1126* 1127* test the input parameters. 1128* 1129 info = 0 1130 if( m.lt.0 ) then 1131 info = -1 1132 else if( n.lt.0 ) then 1133 info = -2 1134 else if( lda.lt.max( 1, m ) ) then 1135 info = -4 1136 end if 1137 if( info.ne.0 ) then 1138 call xerbla( 'dgetrf', -info ) 1139 return 1140 end if 1141* 1142* quick return if possible 1143* 1144 if( m.eq.0 .or. n.eq.0 ) 1145 $ return 1146* 1147* determine the block size for this environment. 1148* 1149 nb = ilaenv( 1, 'dgetrf', ' ', m, n, -1, -1 ) 1150 if( nb.le.1 .or. nb.ge.min( m, n ) ) then 1151* 1152* use unblocked code. 1153* 1154 call dgetf2( m, n, a, lda, ipiv, info ) 1155 else 1156* 1157* use blocked code. 1158* 1159 do 20 j = 1, min( m, n ), nb 1160 jb = min( min( m, n )-j+1, nb ) 1161* 1162* factor diagonal and subdiagonal blocks and test for exact 1163* singularity. 1164* 1165 call dgetf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo ) 1166* 1167* adjust info and the pivot indices. 1168* 1169 if( info.eq.0 .and. iinfo.gt.0 ) 1170 $ info = iinfo + j - 1 1171 do 10 i = j, min( m, j+jb-1 ) 1172 ipiv( i ) = j - 1 + ipiv( i ) 1173 10 continue 1174* 1175* apply interchanges to columns 1:j-1. 1176* 1177 call dlaswp( j-1, a, lda, j, j+jb-1, ipiv, 1 ) 1178* 1179 if( j+jb.le.n ) then 1180* 1181* apply interchanges to columns j+jb:n. 1182* 1183 call dlaswp( n-j-jb+1, a( 1, j+jb ), lda, j, j+jb-1, 1184 $ ipiv, 1 ) 1185* 1186* compute block row of u. 1187* 1188 call dtrsm( 'left', 'lower', 'no transpose', 'unit', jb, 1189 $ n-j-jb+1, one, a( j, j ), lda, a( j, j+jb ), 1190 $ lda ) 1191 if( j+jb.le.m ) then 1192* 1193* update trailing submatrix. 1194* 1195 call dgemm( 'no transpose', 'no transpose', m-j-jb+1, 1196 $ n-j-jb+1, jb, -one, a( j+jb, j ), lda, 1197 $ a( j, j+jb ), lda, one, a( j+jb, j+jb ), 1198 $ lda ) 1199 end if 1200 end if 1201 20 continue 1202 end if 1203 return 1204* 1205* end of dgetrf 1206* 1207 end 1208ccc 1209 subroutine dgetf2( m, n, a, lda, ipiv, info ) 1210* 1211* -- lapack routine (version 3.2) -- 1212* -- lapack is a software package provided by univ. of tennessee, -- 1213* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 1214* november 2006 1215* 1216* .. scalar arguments .. 1217 integer info, lda, m, n 1218* .. 1219* .. array arguments .. 1220 integer ipiv( * ) 1221 double precision a( lda, * ) 1222* .. 1223* 1224* purpose 1225* ======= 1226* 1227* dgetf2 computes an lu factorization of a general m-by-n matrix a 1228* using partial pivoting with row interchanges. 1229* 1230* the factorization has the form 1231* a = p * l * u 1232* where p is a permutation matrix, l is lower triangular with unit 1233* diagonal elements (lower trapezoidal if m > n), and u is upper 1234* triangular (upper trapezoidal if m < n). 1235* 1236* this is the right-looking level 2 blas version of the algorithm. 1237* 1238* arguments 1239* ========= 1240* 1241* m (input) integer 1242* the number of rows of the matrix a. m >= 0. 1243* 1244* n (input) integer 1245* the number of columns of the matrix a. n >= 0. 1246* 1247* a (input/output) double precision array, dimension (lda,n) 1248* on entry, the m by n matrix to be factored. 1249* on exit, the factors l and u from the factorization 1250* a = p*l*u; the unit diagonal elements of l are not stored. 1251* 1252* lda (input) integer 1253* the leading dimension of the array a. lda >= max(1,m). 1254* 1255* ipiv (output) integer array, dimension (min(m,n)) 1256* the pivot indices; for 1 <= i <= min(m,n), row i of the 1257* matrix was interchanged with row ipiv(i). 1258* 1259* info (output) integer 1260* = 0: successful exit 1261* < 0: if info = -k, the k-th argument had an illegal value 1262* > 0: if info = k, u(k,k) is exactly zero. the factorization 1263* has been completed, but the factor u is exactly 1264* singular, and division by zero will occur if it is used 1265* to solve a system of equations. 1266* 1267* ===================================================================== 1268* 1269* .. parameters .. 1270 double precision one, zero 1271 parameter ( one = 1.0d+0, zero = 0.0d+0 ) 1272* .. 1273* .. local scalars .. 1274 double precision sfmin 1275 integer i, j, jp 1276* .. 1277* .. external functions .. 1278 double precision dlamch 1279 integer idamax 1280 external dlamch, idamax 1281* .. 1282* .. external subroutines .. 1283 external dger, dscal, dswap, xerbla 1284* .. 1285* .. intrinsic functions .. 1286 intrinsic max, min 1287* .. 1288* .. executable statements .. 1289* 1290* test the input parameters. 1291* 1292 info = 0 1293 if( m.lt.0 ) then 1294 info = -1 1295 else if( n.lt.0 ) then 1296 info = -2 1297 else if( lda.lt.max( 1, m ) ) then 1298 info = -4 1299 end if 1300 if( info.ne.0 ) then 1301 call xerbla( 'dgetf2', -info ) 1302 return 1303 end if 1304* 1305* quick return if possible 1306* 1307 if( m.eq.0 .or. n.eq.0 ) 1308 $ return 1309* 1310* compute machine safe minimum 1311* 1312 sfmin = dlamch('s') 1313* 1314 do 10 j = 1, min( m, n ) 1315* 1316* find pivot and test for singularity. 1317* 1318 jp = j - 1 + idamax( m-j+1, a( j, j ), 1 ) 1319 ipiv( j ) = jp 1320 if( a( jp, j ).ne.zero ) then 1321* 1322* apply the interchange to columns 1:n. 1323* 1324 if( jp.ne.j ) 1325 $ call dswap( n, a( j, 1 ), lda, a( jp, 1 ), lda ) 1326* 1327* compute elements j+1:m of j-th column. 1328* 1329 if( j.lt.m ) then 1330 if( abs(a( j, j )) .ge. sfmin ) then 1331 call dscal( m-j, one / a( j, j ), a( j+1, j ), 1 ) 1332 else 1333 do 20 i = 1, m-j 1334 a( j+i, j ) = a( j+i, j ) / a( j, j ) 1335 20 continue 1336 end if 1337 end if 1338* 1339 else if( info.eq.0 ) then 1340* 1341 info = j 1342 end if 1343* 1344 if( j.lt.min( m, n ) ) then 1345* 1346* update trailing submatrix. 1347* 1348 call dger( m-j, n-j, -one, a( j+1, j ), 1, a( j, j+1 ), lda, 1349 $ a( j+1, j+1 ), lda ) 1350 end if 1351 10 continue 1352 return 1353* 1354* end of dgetf2 1355* 1356 end 1357ccc 1358 subroutine dger(m,n,alpha,x,incx,y,incy,a,lda) 1359* .. scalar arguments .. 1360 double precision alpha 1361 integer incx,incy,lda,m,n 1362* .. 1363* .. array arguments .. 1364 double precision a(lda,*),x(*),y(*) 1365* .. 1366* 1367* purpose 1368* ======= 1369* 1370* dger performs the rank 1 operation 1371* 1372* a := alpha*x*y' + a, 1373* 1374* where alpha is a scalar, x is an m element vector, y is an n element 1375* vector and a is an m by n matrix. 1376* 1377* arguments 1378* ========== 1379* 1380* m - integer. 1381* on entry, m specifies the number of rows of the matrix a. 1382* m must be at least zero. 1383* unchanged on exit. 1384* 1385* n - integer. 1386* on entry, n specifies the number of columns of the matrix a. 1387* n must be at least zero. 1388* unchanged on exit. 1389* 1390* alpha - double precision. 1391* on entry, alpha specifies the scalar alpha. 1392* unchanged on exit. 1393* 1394* x - double precision array of dimension at least 1395* ( 1 + ( m - 1 )*abs( incx ) ). 1396* before entry, the incremented array x must contain the m 1397* element vector x. 1398* unchanged on exit. 1399* 1400* incx - integer. 1401* on entry, incx specifies the increment for the elements of 1402* x. incx must not be zero. 1403* unchanged on exit. 1404* 1405* y - double precision array of dimension at least 1406* ( 1 + ( n - 1 )*abs( incy ) ). 1407* before entry, the incremented array y must contain the n 1408* element vector y. 1409* unchanged on exit. 1410* 1411* incy - integer. 1412* on entry, incy specifies the increment for the elements of 1413* y. incy must not be zero. 1414* unchanged on exit. 1415* 1416* a - double precision array of dimension ( lda, n ). 1417* before entry, the leading m by n part of the array a must 1418* contain the matrix of coefficients. on exit, a is 1419* overwritten by the updated matrix. 1420* 1421* lda - integer. 1422* on entry, lda specifies the first dimension of a as declared 1423* in the calling (sub) program. lda must be at least 1424* max( 1, m ). 1425* unchanged on exit. 1426* 1427* further details 1428* =============== 1429* 1430* level 2 blas routine. 1431* 1432* -- written on 22-october-1986. 1433* jack dongarra, argonne national lab. 1434* jeremy du croz, nag central office. 1435* sven hammarling, nag central office. 1436* richard hanson, sandia national labs. 1437* 1438* ===================================================================== 1439* 1440* .. parameters .. 1441 double precision zero 1442 parameter (zero=0.0d+0) 1443* .. 1444* .. local scalars .. 1445 double precision temp 1446 integer i,info,ix,j,jy,kx 1447* .. 1448* .. external subroutines .. 1449 external xerbla 1450* .. 1451* .. intrinsic functions .. 1452 intrinsic max 1453* .. 1454* 1455* test the input parameters. 1456* 1457 info = 0 1458 if (m.lt.0) then 1459 info = 1 1460 else if (n.lt.0) then 1461 info = 2 1462 else if (incx.eq.0) then 1463 info = 5 1464 else if (incy.eq.0) then 1465 info = 7 1466 else if (lda.lt.max(1,m)) then 1467 info = 9 1468 end if 1469 if (info.ne.0) then 1470 call xerbla('dger ',info) 1471 return 1472 end if 1473* 1474* quick return if possible. 1475* 1476 if ((m.eq.0) .or. (n.eq.0) .or. (alpha.eq.zero)) return 1477* 1478* start the operations. in this version the elements of a are 1479* accessed sequentially with one pass through a. 1480* 1481 if (incy.gt.0) then 1482 jy = 1 1483 else 1484 jy = 1 - (n-1)*incy 1485 end if 1486 if (incx.eq.1) then 1487 do 20 j = 1,n 1488 if (y(jy).ne.zero) then 1489 temp = alpha*y(jy) 1490 do 10 i = 1,m 1491 a(i,j) = a(i,j) + x(i)*temp 1492 10 continue 1493 end if 1494 jy = jy + incy 1495 20 continue 1496 else 1497 if (incx.gt.0) then 1498 kx = 1 1499 else 1500 kx = 1 - (m-1)*incx 1501 end if 1502 do 40 j = 1,n 1503 if (y(jy).ne.zero) then 1504 temp = alpha*y(jy) 1505 ix = kx 1506 do 30 i = 1,m 1507 a(i,j) = a(i,j) + x(ix)*temp 1508 ix = ix + incx 1509 30 continue 1510 end if 1511 jy = jy + incy 1512 40 continue 1513 end if 1514* 1515 return 1516* 1517* end of dger . 1518* 1519 end 1520ccc 1521 subroutine dscal(n,da,dx,incx) 1522* .. scalar arguments .. 1523 double precision da 1524 integer incx,n 1525* .. 1526* .. array arguments .. 1527 double precision dx(*) 1528* .. 1529* 1530* purpose 1531* ======= 1532* 1533* dscal scales a vector by a constant. 1534* uses unrolled loops for increment equal to one. 1535* 1536* further details 1537* =============== 1538* 1539* jack dongarra, linpack, 3/11/78. 1540* modified 3/93 to return if incx .le. 0. 1541* modified 12/3/93, array(1) declarations changed to array(*) 1542* 1543* ===================================================================== 1544* 1545* .. local scalars .. 1546 integer i,m,mp1,nincx 1547* .. 1548* .. intrinsic functions .. 1549 intrinsic mod 1550* .. 1551 if (n.le.0 .or. incx.le.0) return 1552 if (incx.eq.1) go to 20 1553* 1554* code for increment not equal to 1 1555* 1556 nincx = n*incx 1557 do 10 i = 1,nincx,incx 1558 dx(i) = da*dx(i) 1559 10 continue 1560 return 1561* 1562* code for increment equal to 1 1563* 1564* 1565* clean-up loop 1566* 1567 20 m = mod(n,5) 1568 if (m.eq.0) go to 40 1569 do 30 i = 1,m 1570 dx(i) = da*dx(i) 1571 30 continue 1572 if (n.lt.5) return 1573 40 mp1 = m + 1 1574 do 50 i = mp1,n,5 1575 dx(i) = da*dx(i) 1576 dx(i+1) = da*dx(i+1) 1577 dx(i+2) = da*dx(i+2) 1578 dx(i+3) = da*dx(i+3) 1579 dx(i+4) = da*dx(i+4) 1580 50 continue 1581 return 1582 end 1583ccc 1584 subroutine dswap(n,dx,incx,dy,incy) 1585* .. scalar arguments .. 1586 integer incx,incy,n 1587* .. 1588* .. array arguments .. 1589 double precision dx(*),dy(*) 1590* .. 1591* 1592* purpose 1593* ======= 1594* 1595* interchanges two vectors. 1596* uses unrolled loops for increments equal one. 1597* 1598* further details 1599* =============== 1600* 1601* jack dongarra, linpack, 3/11/78. 1602* modified 12/3/93, array(1) declarations changed to array(*) 1603* 1604* ===================================================================== 1605* 1606* .. local scalars .. 1607 double precision dtemp 1608 integer i,ix,iy,m,mp1 1609* .. 1610* .. intrinsic functions .. 1611 intrinsic mod 1612* .. 1613 if (n.le.0) return 1614 if (incx.eq.1 .and. incy.eq.1) go to 20 1615* 1616* code for unequal increments or equal increments not equal 1617* to 1 1618* 1619 ix = 1 1620 iy = 1 1621 if (incx.lt.0) ix = (-n+1)*incx + 1 1622 if (incy.lt.0) iy = (-n+1)*incy + 1 1623 do 10 i = 1,n 1624 dtemp = dx(ix) 1625 dx(ix) = dy(iy) 1626 dy(iy) = dtemp 1627 ix = ix + incx 1628 iy = iy + incy 1629 10 continue 1630 return 1631* 1632* code for both increments equal to 1 1633* 1634* 1635* clean-up loop 1636* 1637 20 m = mod(n,3) 1638 if (m.eq.0) go to 40 1639 do 30 i = 1,m 1640 dtemp = dx(i) 1641 dx(i) = dy(i) 1642 dy(i) = dtemp 1643 30 continue 1644 if (n.lt.3) return 1645 40 mp1 = m + 1 1646 do 50 i = mp1,n,3 1647 dtemp = dx(i) 1648 dx(i) = dy(i) 1649 dy(i) = dtemp 1650 dtemp = dx(i+1) 1651 dx(i+1) = dy(i+1) 1652 dy(i+1) = dtemp 1653 dtemp = dx(i+2) 1654 dx(i+2) = dy(i+2) 1655 dy(i+2) = dtemp 1656 50 continue 1657 return 1658 end 1659ccc 1660 subroutine xerbla( srname, info ) 1661* 1662* -- lapack auxiliary routine (preliminary version) -- 1663* univ. of tennessee, univ. of california berkeley and nag ltd.. 1664* november 2006 1665* 1666* .. scalar arguments .. 1667 character*(*) srname 1668 integer info 1669* .. 1670* 1671* purpose 1672* ======= 1673* 1674* xerbla is an error handler for the lapack routines. 1675* it is called by an lapack routine if an input parameter has an 1676* invalid value. a message is printed and execution stops. 1677* 1678* installers may consider modifying the stop statement in order to 1679* call system-specific exception-handling facilities. 1680* 1681* arguments 1682* ========= 1683* 1684* srname (input) character*(*) 1685* the name of the routine which called xerbla. 1686* 1687* info (input) integer 1688* the position of the invalid parameter in the parameter list 1689* of the calling routine. 1690* 1691* ===================================================================== 1692* 1693* .. intrinsic functions .. 1694 intrinsic len_trim 1695* .. 1696* .. executable statements .. 1697* 1698 write( *, fmt = 9999 )srname( 1:len_trim( srname ) ), info 1699* 1700 stop 1701* 1702 9999 format( ' ** on entry to ', a, ' parameter number ', i2, ' had ', 1703 $ 'an illegal value' ) 1704* 1705* end of xerbla 1706* 1707 end 1708ccc 1709 subroutine dgetrs( trans, n, nrhs, a, lda, ipiv, b, ldb, info ) 1710* 1711* -- lapack routine (version 3.2) -- 1712* -- lapack is a software package provided by univ. of tennessee, -- 1713* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 1714* november 2006 1715* 1716* .. scalar arguments .. 1717 character trans 1718 integer info, lda, ldb, n, nrhs 1719* .. 1720* .. array arguments .. 1721 integer ipiv( * ) 1722 double precision a( lda, * ), b( ldb, * ) 1723* .. 1724* 1725* purpose 1726* ======= 1727* 1728* dgetrs solves a system of linear equations 1729* a * x = b or a' * x = b 1730* with a general n-by-n matrix a using the lu factorization computed 1731* by dgetrf. 1732* 1733* arguments 1734* ========= 1735* 1736* trans (input) character*1 1737* specifies the form of the system of equations: 1738* = 'n': a * x = b (no transpose) 1739* = 't': a'* x = b (transpose) 1740* = 'c': a'* x = b (conjugate transpose = transpose) 1741* 1742* n (input) integer 1743* the order of the matrix a. n >= 0. 1744* 1745* nrhs (input) integer 1746* the number of right hand sides, i.e., the number of columns 1747* of the matrix b. nrhs >= 0. 1748* 1749* a (input) double precision array, dimension (lda,n) 1750* the factors l and u from the factorization a = p*l*u 1751* as computed by dgetrf. 1752* 1753* lda (input) integer 1754* the leading dimension of the array a. lda >= max(1,n). 1755* 1756* ipiv (input) integer array, dimension (n) 1757* the pivot indices from dgetrf; for 1<=i<=n, row i of the 1758* matrix was interchanged with row ipiv(i). 1759* 1760* b (input/output) double precision array, dimension (ldb,nrhs) 1761* on entry, the right hand side matrix b. 1762* on exit, the solution matrix x. 1763* 1764* ldb (input) integer 1765* the leading dimension of the array b. ldb >= max(1,n). 1766* 1767* info (output) integer 1768* = 0: successful exit 1769* < 0: if info = -i, the i-th argument had an illegal value 1770* 1771* ===================================================================== 1772* 1773* .. parameters .. 1774 double precision one 1775 parameter ( one = 1.0d+0 ) 1776* .. 1777* .. local scalars .. 1778 logical notran 1779* .. 1780* .. external functions .. 1781 logical lsame 1782 external lsame 1783* .. 1784* .. external subroutines .. 1785 external dlaswp, dtrsm, xerbla 1786* .. 1787* .. intrinsic functions .. 1788 intrinsic max 1789* .. 1790* .. executable statements .. 1791* 1792* test the input parameters. 1793* 1794 info = 0 1795 notran = lsame( trans, 'n' ) 1796 if( .not.notran .and. .not.lsame( trans, 't' ) .and. .not. 1797 $ lsame( trans, 'c' ) ) then 1798 info = -1 1799 else if( n.lt.0 ) then 1800 info = -2 1801 else if( nrhs.lt.0 ) then 1802 info = -3 1803 else if( lda.lt.max( 1, n ) ) then 1804 info = -5 1805 else if( ldb.lt.max( 1, n ) ) then 1806 info = -8 1807 end if 1808 if( info.ne.0 ) then 1809 call xerbla( 'dgetrs', -info ) 1810 return 1811 end if 1812* 1813* quick return if possible 1814* 1815 if( n.eq.0 .or. nrhs.eq.0 ) 1816 $ return 1817* 1818 if( notran ) then 1819* 1820* solve a * x = b. 1821* 1822* apply row interchanges to the right hand sides. 1823* 1824 call dlaswp( nrhs, b, ldb, 1, n, ipiv, 1 ) 1825* 1826* solve l*x = b, overwriting b with x. 1827* 1828 call dtrsm( 'left', 'lower', 'no transpose', 'unit', n, nrhs, 1829 $ one, a, lda, b, ldb ) 1830* 1831* solve u*x = b, overwriting b with x. 1832* 1833 call dtrsm( 'left', 'upper', 'no transpose', 'non-unit', n, 1834 $ nrhs, one, a, lda, b, ldb ) 1835 else 1836* 1837* solve a' * x = b. 1838* 1839* solve u'*x = b, overwriting b with x. 1840* 1841 call dtrsm( 'left', 'upper', 'transpose', 'non-unit', n, nrhs, 1842 $ one, a, lda, b, ldb ) 1843* 1844* solve l'*x = b, overwriting b with x. 1845* 1846 call dtrsm( 'left', 'lower', 'transpose', 'unit', n, nrhs, one, 1847 $ a, lda, b, ldb ) 1848* 1849* apply row interchanges to the solution vectors. 1850* 1851 call dlaswp( nrhs, b, ldb, 1, n, ipiv, -1 ) 1852 end if 1853* 1854 return 1855* 1856* end of dgetrs 1857* 1858 end 1859ccc 1860 subroutine dgemm(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc) 1861* .. scalar arguments .. 1862 double precision alpha,beta 1863 integer k,lda,ldb,ldc,m,n 1864 character transa,transb 1865* .. 1866* .. array arguments .. 1867 double precision a(lda,*),b(ldb,*),c(ldc,*) 1868* .. 1869* 1870* purpose 1871* ======= 1872* 1873* dgemm performs one of the matrix-matrix operations 1874* 1875* c := alpha*op( a )*op( b ) + beta*c, 1876* 1877* where op( x ) is one of 1878* 1879* op( x ) = x or op( x ) = x', 1880* 1881* alpha and beta are scalars, and a, b and c are matrices, with op( a ) 1882* an m by k matrix, op( b ) a k by n matrix and c an m by n matrix. 1883* 1884* arguments 1885* ========== 1886* 1887* transa - character*1. 1888* on entry, transa specifies the form of op( a ) to be used in 1889* the matrix multiplication as follows: 1890* 1891* transa = 'n' or 'n', op( a ) = a. 1892* 1893* transa = 't' or 't', op( a ) = a'. 1894* 1895* transa = 'c' or 'c', op( a ) = a'. 1896* 1897* unchanged on exit. 1898* 1899* transb - character*1. 1900* on entry, transb specifies the form of op( b ) to be used in 1901* the matrix multiplication as follows: 1902* 1903* transb = 'n' or 'n', op( b ) = b. 1904* 1905* transb = 't' or 't', op( b ) = b'. 1906* 1907* transb = 'c' or 'c', op( b ) = b'. 1908* 1909* unchanged on exit. 1910* 1911* m - integer. 1912* on entry, m specifies the number of rows of the matrix 1913* op( a ) and of the matrix c. m must be at least zero. 1914* unchanged on exit. 1915* 1916* n - integer. 1917* on entry, n specifies the number of columns of the matrix 1918* op( b ) and the number of columns of the matrix c. n must be 1919* at least zero. 1920* unchanged on exit. 1921* 1922* k - integer. 1923* on entry, k specifies the number of columns of the matrix 1924* op( a ) and the number of rows of the matrix op( b ). k must 1925* be at least zero. 1926* unchanged on exit. 1927* 1928* alpha - double precision. 1929* on entry, alpha specifies the scalar alpha. 1930* unchanged on exit. 1931* 1932* a - double precision array of dimension ( lda, ka ), where ka is 1933* k when transa = 'n' or 'n', and is m otherwise. 1934* before entry with transa = 'n' or 'n', the leading m by k 1935* part of the array a must contain the matrix a, otherwise 1936* the leading k by m part of the array a must contain the 1937* matrix a. 1938* unchanged on exit. 1939* 1940* lda - integer. 1941* on entry, lda specifies the first dimension of a as declared 1942* in the calling (sub) program. when transa = 'n' or 'n' then 1943* lda must be at least max( 1, m ), otherwise lda must be at 1944* least max( 1, k ). 1945* unchanged on exit. 1946* 1947* b - double precision array of dimension ( ldb, kb ), where kb is 1948* n when transb = 'n' or 'n', and is k otherwise. 1949* before entry with transb = 'n' or 'n', the leading k by n 1950* part of the array b must contain the matrix b, otherwise 1951* the leading n by k part of the array b must contain the 1952* matrix b. 1953* unchanged on exit. 1954* 1955* ldb - integer. 1956* on entry, ldb specifies the first dimension of b as declared 1957* in the calling (sub) program. when transb = 'n' or 'n' then 1958* ldb must be at least max( 1, k ), otherwise ldb must be at 1959* least max( 1, n ). 1960* unchanged on exit. 1961* 1962* beta - double precision. 1963* on entry, beta specifies the scalar beta. when beta is 1964* supplied as zero then c need not be set on input. 1965* unchanged on exit. 1966* 1967* c - double precision array of dimension ( ldc, n ). 1968* before entry, the leading m by n part of the array c must 1969* contain the matrix c, except when beta is zero, in which 1970* case c need not be set on entry. 1971* on exit, the array c is overwritten by the m by n matrix 1972* ( alpha*op( a )*op( b ) + beta*c ). 1973* 1974* ldc - integer. 1975* on entry, ldc specifies the first dimension of c as declared 1976* in the calling (sub) program. ldc must be at least 1977* max( 1, m ). 1978* unchanged on exit. 1979* 1980* further details 1981* =============== 1982* 1983* level 3 blas routine. 1984* 1985* -- written on 8-february-1989. 1986* jack dongarra, argonne national laboratory. 1987* iain duff, aere harwell. 1988* jeremy du croz, numerical algorithms group ltd. 1989* sven hammarling, numerical algorithms group ltd. 1990* 1991* ===================================================================== 1992* 1993* .. external functions .. 1994 logical lsame 1995 external lsame 1996* .. 1997* .. external subroutines .. 1998 external xerbla 1999* .. 2000* .. intrinsic functions .. 2001 intrinsic max 2002* .. 2003* .. local scalars .. 2004 double precision temp 2005 integer i,info,j,l,ncola,nrowa,nrowb 2006 logical nota,notb 2007* .. 2008* .. parameters .. 2009 double precision one,zero 2010 parameter (one=1.0d+0,zero=0.0d+0) 2011* .. 2012* 2013* set nota and notb as true if a and b respectively are not 2014* transposed and set nrowa, ncola and nrowb as the number of rows 2015* and columns of a and the number of rows of b respectively. 2016* 2017 nota = lsame(transa,'n') 2018 notb = lsame(transb,'n') 2019 if (nota) then 2020 nrowa = m 2021 ncola = k 2022 else 2023 nrowa = k 2024 ncola = m 2025 end if 2026 if (notb) then 2027 nrowb = k 2028 else 2029 nrowb = n 2030 end if 2031* 2032* test the input parameters. 2033* 2034 info = 0 2035 if ((.not.nota) .and. (.not.lsame(transa,'c')) .and. 2036 + (.not.lsame(transa,'t'))) then 2037 info = 1 2038 else if ((.not.notb) .and. (.not.lsame(transb,'c')) .and. 2039 + (.not.lsame(transb,'t'))) then 2040 info = 2 2041 else if (m.lt.0) then 2042 info = 3 2043 else if (n.lt.0) then 2044 info = 4 2045 else if (k.lt.0) then 2046 info = 5 2047 else if (lda.lt.max(1,nrowa)) then 2048 info = 8 2049 else if (ldb.lt.max(1,nrowb)) then 2050 info = 10 2051 else if (ldc.lt.max(1,m)) then 2052 info = 13 2053 end if 2054 if (info.ne.0) then 2055 call xerbla('dgemm ',info) 2056 return 2057 end if 2058* 2059* quick return if possible. 2060* 2061 if ((m.eq.0) .or. (n.eq.0) .or. 2062 + (((alpha.eq.zero).or. (k.eq.0)).and. (beta.eq.one))) return 2063* 2064* and if alpha.eq.zero. 2065* 2066 if (alpha.eq.zero) then 2067 if (beta.eq.zero) then 2068 do 20 j = 1,n 2069 do 10 i = 1,m 2070 c(i,j) = zero 2071 10 continue 2072 20 continue 2073 else 2074 do 40 j = 1,n 2075 do 30 i = 1,m 2076 c(i,j) = beta*c(i,j) 2077 30 continue 2078 40 continue 2079 end if 2080 return 2081 end if 2082* 2083* start the operations. 2084* 2085 if (notb) then 2086 if (nota) then 2087* 2088* form c := alpha*a*b + beta*c. 2089* 2090 do 90 j = 1,n 2091 if (beta.eq.zero) then 2092 do 50 i = 1,m 2093 c(i,j) = zero 2094 50 continue 2095 else if (beta.ne.one) then 2096 do 60 i = 1,m 2097 c(i,j) = beta*c(i,j) 2098 60 continue 2099 end if 2100 do 80 l = 1,k 2101 if (b(l,j).ne.zero) then 2102 temp = alpha*b(l,j) 2103 do 70 i = 1,m 2104 c(i,j) = c(i,j) + temp*a(i,l) 2105 70 continue 2106 end if 2107 80 continue 2108 90 continue 2109 else 2110* 2111* form c := alpha*a'*b + beta*c 2112* 2113 do 120 j = 1,n 2114 do 110 i = 1,m 2115 temp = zero 2116 do 100 l = 1,k 2117 temp = temp + a(l,i)*b(l,j) 2118 100 continue 2119 if (beta.eq.zero) then 2120 c(i,j) = alpha*temp 2121 else 2122 c(i,j) = alpha*temp + beta*c(i,j) 2123 end if 2124 110 continue 2125 120 continue 2126 end if 2127 else 2128 if (nota) then 2129* 2130* form c := alpha*a*b' + beta*c 2131* 2132 do 170 j = 1,n 2133 if (beta.eq.zero) then 2134 do 130 i = 1,m 2135 c(i,j) = zero 2136 130 continue 2137 else if (beta.ne.one) then 2138 do 140 i = 1,m 2139 c(i,j) = beta*c(i,j) 2140 140 continue 2141 end if 2142 do 160 l = 1,k 2143 if (b(j,l).ne.zero) then 2144 temp = alpha*b(j,l) 2145 do 150 i = 1,m 2146 c(i,j) = c(i,j) + temp*a(i,l) 2147 150 continue 2148 end if 2149 160 continue 2150 170 continue 2151 else 2152* 2153* form c := alpha*a'*b' + beta*c 2154* 2155 do 200 j = 1,n 2156 do 190 i = 1,m 2157 temp = zero 2158 do 180 l = 1,k 2159 temp = temp + a(l,i)*b(j,l) 2160 180 continue 2161 if (beta.eq.zero) then 2162 c(i,j) = alpha*temp 2163 else 2164 c(i,j) = alpha*temp + beta*c(i,j) 2165 end if 2166 190 continue 2167 200 continue 2168 end if 2169 end if 2170* 2171 return 2172* 2173* end of dgemm . 2174* 2175 end 2176ccc 2177 subroutine dtrsm(side,uplo,transa,diag,m,n,alpha,a,lda,b,ldb) 2178* .. scalar arguments .. 2179 double precision alpha 2180 integer lda,ldb,m,n 2181 character diag,side,transa,uplo 2182* .. 2183* .. array arguments .. 2184 double precision a(lda,*),b(ldb,*) 2185* .. 2186* 2187* purpose 2188* ======= 2189* 2190* dtrsm solves one of the matrix equations 2191* 2192* op( a )*x = alpha*b, or x*op( a ) = alpha*b, 2193* 2194* where alpha is a scalar, x and b are m by n matrices, a is a unit, or 2195* non-unit, upper or lower triangular matrix and op( a ) is one of 2196* 2197* op( a ) = a or op( a ) = a'. 2198* 2199* the matrix x is overwritten on b. 2200* 2201* arguments 2202* ========== 2203* 2204* side - character*1. 2205* on entry, side specifies whether op( a ) appears on the left 2206* or right of x as follows: 2207* 2208* side = 'l' or 'l' op( a )*x = alpha*b. 2209* 2210* side = 'r' or 'r' x*op( a ) = alpha*b. 2211* 2212* unchanged on exit. 2213* 2214* uplo - character*1. 2215* on entry, uplo specifies whether the matrix a is an upper or 2216* lower triangular matrix as follows: 2217* 2218* uplo = 'u' or 'u' a is an upper triangular matrix. 2219* 2220* uplo = 'l' or 'l' a is a lower triangular matrix. 2221* 2222* unchanged on exit. 2223* 2224* transa - character*1. 2225* on entry, transa specifies the form of op( a ) to be used in 2226* the matrix multiplication as follows: 2227* 2228* transa = 'n' or 'n' op( a ) = a. 2229* 2230* transa = 't' or 't' op( a ) = a'. 2231* 2232* transa = 'c' or 'c' op( a ) = a'. 2233* 2234* unchanged on exit. 2235* 2236* diag - character*1. 2237* on entry, diag specifies whether or not a is unit triangular 2238* as follows: 2239* 2240* diag = 'u' or 'u' a is assumed to be unit triangular. 2241* 2242* diag = 'n' or 'n' a is not assumed to be unit 2243* triangular. 2244* 2245* unchanged on exit. 2246* 2247* m - integer. 2248* on entry, m specifies the number of rows of b. m must be at 2249* least zero. 2250* unchanged on exit.gchuev-874 2251* 2252* n - integer. 2253* on entry, n specifies the number of columns of b. n must be 2254* at least zero. 2255* unchanged on exit. 2256* 2257* alpha - double precision. 2258* on entry, alpha specifies the scalar alpha. when alpha is 2259* zero then a is not referenced and b need not be set before 2260* entry. 2261* unchanged on exit. 2262* 2263* a - double precision array of dimension ( lda, k ), where k is m 2264* when side = 'l' or 'l' and is n when side = 'r' or 'r'. 2265* before entry with uplo = 'u' or 'u', the leading k by k 2266* upper triangular part of the array a must contain the upper 2267* triangular matrix and the strictly lower triangular part of 2268* a is not referenced. 2269* before entry with uplo = 'l' or 'l', the leading k by k 2270* lower triangular part of the array a must contain the lower 2271* triangular matrix and the strictly upper triangular part of 2272* a is not referenced. 2273* note that when diag = 'u' or 'u', the diagonal elements of 2274* a are not referenced either, but are assumed to be unity. 2275* unchanged on exit. 2276* 2277* lda - integer. 2278* on entry, lda specifies the first dimension of a as declared 2279* in the calling (sub) program. when side = 'l' or 'l' then 2280* lda must be at least max( 1, m ), when side = 'r' or 'r' 2281* then lda must be at least max( 1, n ). 2282* unchanged on exit. 2283* 2284* b - double precision array of dimension ( ldb, n ). 2285* before entry, the leading m by n part of the array b must 2286* contain the right-hand side matrix b, and on exit is 2287* overwritten by the solution matrix x. 2288* 2289* ldb - integer. 2290* on entry, ldb specifies the first dimension of b as declared 2291* in the calling (sub) program. ldb must be at least 2292* max( 1, m ). 2293* unchanged on exit. 2294* 2295* further details 2296* =============== 2297* 2298* level 3 blas routine. 2299* 2300* 2301* -- written on 8-february-1989. 2302* jack dongarra, argonne national laboratory. 2303* iain duff, aere harwell. 2304* jeremy du croz, numerical algorithms group ltd. 2305* sven hammarling, numerical algorithms group ltd. 2306* 2307* ===================================================================== 2308* 2309* .. external functions .. 2310 logical lsame 2311 external lsame 2312* .. 2313* .. external subroutines .. 2314 external xerbla 2315* .. 2316* .. intrinsic functions .. 2317 intrinsic max 2318* .. 2319* .. local scalars .. 2320 double precision temp 2321 integer i,info,j,k,nrowa 2322 logical lside,nounit,upper 2323* .. 2324* .. parameters .. 2325 double precision one,zero 2326 parameter (one=1.0d+0,zero=0.0d+0) 2327* .. 2328* 2329* test the input parameters. 2330* 2331 lside = lsame(side,'l') 2332 if (lside) then 2333 nrowa = m 2334 else 2335 nrowa = n 2336 end if 2337 nounit = lsame(diag,'n') 2338 upper = lsame(uplo,'u') 2339* 2340 info = 0 2341 if ((.not.lside) .and. (.not.lsame(side,'r'))) then 2342 info = 1 2343 else if ((.not.upper) .and. (.not.lsame(uplo,'l'))) then 2344 info = 2 2345 else if ((.not.lsame(transa,'n')) .and. 2346 + (.not.lsame(transa,'t')) .and. 2347 + (.not.lsame(transa,'c'))) then 2348 info = 3 2349 else if ((.not.lsame(diag,'u')) .and. (.not.lsame(diag,'n'))) then 2350 info = 4 2351 else if (m.lt.0) then 2352 info = 5 2353 else if (n.lt.0) then 2354 info = 6 2355 else if (lda.lt.max(1,nrowa)) then 2356 info = 9 2357 else if (ldb.lt.max(1,m)) then 2358 info = 11 2359 end if 2360 if (info.ne.0) then 2361 call xerbla('dtrsm ',info) 2362 return 2363 end if 2364* 2365* quick return if possible. 2366* 2367 if (m.eq.0 .or. n.eq.0) return 2368* 2369* and when alpha.eq.zero. 2370* 2371 if (alpha.eq.zero) then 2372 do 20 j = 1,n 2373 do 10 i = 1,m 2374 b(i,j) = zero 2375 10 continue 2376 20 continue 2377 return 2378 end if 2379* 2380* start the operations. 2381* 2382 if (lside) then 2383 if (lsame(transa,'n')) then 2384* 2385* form b := alpha*inv( a )*b. 2386* 2387 if (upper) then 2388 do 60 j = 1,n 2389 if (alpha.ne.one) then 2390 do 30 i = 1,m 2391 b(i,j) = alpha*b(i,j) 2392 30 continue 2393 end if 2394 do 50 k = m,1,-1 2395 if (b(k,j).ne.zero) then 2396 if (nounit) b(k,j) = b(k,j)/a(k,k) 2397 do 40 i = 1,k - 1 2398 b(i,j) = b(i,j) - b(k,j)*a(i,k) 2399 40 continue 2400 end if 2401 50 continue 2402 60 continue 2403 else 2404 do 100 j = 1,n 2405 if (alpha.ne.one) then 2406 do 70 i = 1,m 2407 b(i,j) = alpha*b(i,j) 2408 70 continue 2409 end if 2410 do 90 k = 1,m 2411 if (b(k,j).ne.zero) then 2412 if (nounit) b(k,j) = b(k,j)/a(k,k) 2413 do 80 i = k + 1,m 2414 b(i,j) = b(i,j) - b(k,j)*a(i,k) 2415 80 continue 2416 end if 2417 90 continue 2418 100 continue 2419 end if 2420 else 2421* 2422* form b := alpha*inv( a' )*b. 2423* 2424 if (upper) then 2425 do 130 j = 1,n 2426 do 120 i = 1,m 2427 temp = alpha*b(i,j) 2428 do 110 k = 1,i - 1 2429 temp = temp - a(k,i)*b(k,j) 2430 110 continue 2431 if (nounit) temp = temp/a(i,i) 2432 b(i,j) = temp 2433 120 continue 2434 130 continue 2435 else 2436 do 160 j = 1,n 2437 do 150 i = m,1,-1 2438 temp = alpha*b(i,j) 2439 do 140 k = i + 1,m 2440 temp = temp - a(k,i)*b(k,j) 2441 140 continue 2442 if (nounit) temp = temp/a(i,i) 2443 b(i,j) = temp 2444 150 continue 2445 160 continue 2446 end if 2447 end if 2448 else 2449 if (lsame(transa,'n')) then 2450* 2451* form b := alpha*b*inv( a ). 2452* 2453 if (upper) then 2454 do 210 j = 1,n 2455 if (alpha.ne.one) then 2456 do 170 i = 1,m 2457 b(i,j) = alpha*b(i,j) 2458 170 continue 2459 end if 2460 do 190 k = 1,j - 1 2461 if (a(k,j).ne.zero) then 2462 do 180 i = 1,m 2463 b(i,j) = b(i,j) - a(k,j)*b(i,k) 2464 180 continue 2465 end if 2466 190 continue 2467 if (nounit) then 2468 temp = one/a(j,j) 2469 do 200 i = 1,m 2470 b(i,j) = temp*b(i,j) 2471 200 continue 2472 end if 2473 210 continue 2474 else 2475 do 260 j = n,1,-1 2476 if (alpha.ne.one) then 2477 do 220 i = 1,m 2478 b(i,j) = alpha*b(i,j) 2479 220 continue 2480 end if 2481 do 240 k = j + 1,n 2482 if (a(k,j).ne.zero) then 2483 do 230 i = 1,m 2484 b(i,j) = b(i,j) - a(k,j)*b(i,k) 2485 230 continue 2486 end if 2487 240 continue 2488 if (nounit) then 2489 temp = one/a(j,j) 2490 do 250 i = 1,m 2491 b(i,j) = temp*b(i,j) 2492 250 continue 2493 end if 2494 260 continue 2495 end if 2496 else 2497* 2498* form b := alpha*b*inv( a' ). 2499* 2500 if (upper) then 2501 do 310 k = n,1,-1 2502 if (nounit) then 2503 temp = one/a(k,k) 2504 do 270 i = 1,m 2505 b(i,k) = temp*b(i,k) 2506 270 continue 2507 end if 2508 do 290 j = 1,k - 1 2509 if (a(j,k).ne.zero) then 2510 temp = a(j,k) 2511 do 280 i = 1,m 2512 b(i,j) = b(i,j) - temp*b(i,k) 2513 280 continue 2514 end if 2515 290 continue 2516 if (alpha.ne.one) then 2517 do 300 i = 1,m 2518 b(i,k) = alpha*b(i,k) 2519 300 continue 2520 end if 2521 310 continue 2522 else 2523 do 360 k = 1,n 2524 if (nounit) then 2525 temp = one/a(k,k) 2526 do 320 i = 1,m 2527 b(i,k) = temp*b(i,k) 2528 320 continue 2529 end if 2530 do 340 j = k + 1,n 2531 if (a(j,k).ne.zero) then 2532 temp = a(j,k) 2533 do 330 i = 1,m 2534 b(i,j) = b(i,j) - temp*b(i,k) 2535 330 continue 2536 end if 2537 340 continue 2538 if (alpha.ne.one) then 2539 do 350 i = 1,m 2540 b(i,k) = alpha*b(i,k) 2541 350 continue 2542 end if 2543 360 continue 2544 end if 2545 end if 2546 end if 2547* 2548 return 2549* 2550* end of dtrsm . 2551* 2552 end 2553ccc 2554 subroutine dlaswp( n, a, lda, k1, k2, ipiv, incx ) 2555* 2556* -- lapack auxiliary routine (version 3.2) -- 2557* -- lapack is a software package provided by univ. of tennessee, -- 2558* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 2559* november 2006 2560* 2561* .. scalar arguments .. 2562 integer incx, k1, k2, lda, n 2563* .. 2564* .. array arguments .. 2565 integer ipiv( * ) 2566 double precision a( lda, * ) 2567* .. 2568* 2569* purpose 2570* ======= 2571* 2572* dlaswp performs a series of row interchanges on the matrix a. 2573* one row interchange is initiated for each of rows k1 through k2 of a. 2574* 2575* arguments 2576* ========= 2577* 2578* n (input) integer 2579* the number of columns of the matrix a. 2580* 2581* a (input/output) double precision array, dimension (lda,n) 2582* on entry, the matrix of column dimension n to which the row 2583* interchanges will be applied. 2584* on exit, the permuted matrix. 2585* 2586* lda (input) integer 2587* the leading dimension of the array a. 2588* 2589* k1 (input) integer 2590* the first element of ipiv for which a row interchange will 2591* be done. 2592* 2593* k2 (input) integer 2594* the last element of ipiv for which a row interchange will 2595* be done. 2596* 2597* ipiv (input) integer array, dimension (k2*abs(incx)) 2598* the vector of pivot indices. only the elements in positions 2599* k1 through k2 of ipiv are accessed. 2600* ipiv(k) = l implies rows k and l are to be interchanged. 2601* 2602* incx (input) integer 2603* the increment between successive values of ipiv. if ipiv 2604* is negative, the pivots are applied in reverse order. 2605* 2606* further details 2607* =============== 2608* 2609* modified by 2610* r. c. whaley, computer science dept., univ. of tenn., knoxville, usa 2611* 2612* ===================================================================== 2613* 2614* .. local scalars .. 2615 integer i, i1, i2, inc, ip, ix, ix0, j, k, n32 2616 double precision temp 2617* .. 2618* .. executable statements .. 2619* 2620* interchange row i with row ipiv(i) for each of rows k1 through k2. 2621* 2622 if( incx.gt.0 ) then 2623 ix0 = k1 2624 i1 = k1 2625 i2 = k2 2626 inc = 1 2627 else if( incx.lt.0 ) then 2628 ix0 = 1 + ( 1-k2 )*incx 2629 i1 = k2 2630 i2 = k1 2631 inc = -1 2632 else 2633 return 2634 end if 2635* 2636 n32 = ( n / 32 )*32 2637 if( n32.ne.0 ) then 2638 do 30 j = 1, n32, 32 2639 ix = ix0 2640 do 20 i = i1, i2, inc 2641 ip = ipiv( ix ) 2642 if( ip.ne.i ) then 2643 do 10 k = j, j + 31 2644 temp = a( i, k ) 2645 a( i, k ) = a( ip, k ) 2646 a( ip, k ) = temp 2647 10 continue 2648 end if 2649 ix = ix + incx 2650 20 continue 2651 30 continue 2652 end if 2653 if( n32.ne.n ) then 2654 n32 = n32 + 1 2655 ix = ix0 2656 do 50 i = i1, i2, inc 2657 ip = ipiv( ix ) 2658 if( ip.ne.i ) then 2659 do 40 k = n32, n 2660 temp = a( i, k ) 2661 a( i, k ) = a( ip, k ) 2662 a( ip, k ) = temp 2663 40 continue 2664 end if 2665 ix = ix + incx 2666 50 continue 2667 end if 2668* 2669 return 2670* 2671* end of dlaswp 2672* 2673 end 2674ccc 2675 integer function ilaenv( ispec, name, opts, n1, n2, n3, n4 ) 2676* 2677* -- lapack auxiliary routine (version 3.2.1) -- 2678* 2679* -- april 2009 -- 2680* 2681* -- lapack is a software package provided by univ. of tennessee, -- 2682* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 2683* 2684* .. scalar arguments .. 2685 character*( * ) name, opts 2686 integer ispec, n1, n2, n3, n4 2687* .. 2688* 2689* purpose 2690* ======= 2691* 2692* ilaenv is called from the lapack routines to choose problem-dependent 2693* parameters for the local environment. see ispec for a description of 2694* the parameters. 2695* 2696* ilaenv returns an integer 2697* if ilaenv >= 0: ilaenv returns the value of the parameter specified by ispec 2698* if ilaenv < 0: if ilaenv = -k, the k-th argument had an illegal value. 2699* 2700* this version provides a set of parameters which should give good, 2701* but not optimal, performance on many of the currently available 2702* computers. users are encouraged to modify this subroutine to set 2703* the tuning parameters for their particular machine using the option 2704* and problem size information in the arguments. 2705* 2706* this routine will not function correctly if it is converted to all 2707* lower case. converting it to all upper case is allowed. 2708* 2709* arguments 2710* ========= 2711* 2712* ispec (input) integer 2713* specifies the parameter to be returned as the value of 2714* ilaenv. 2715* = 1: the optimal blocksize; if this value is 1, an unblocked 2716* algorithm will give the best performance. 2717* = 2: the minimum block size for which the block routine 2718* should be used; if the usable block size is less than 2719* this value, an unblocked routine should be used. 2720* = 3: the crossover point (in a block routine, for n less 2721* than this value, an unblocked routine should be used) 2722* = 4: the number of shifts, used in the nonsymmetric 2723* eigenvalue routines (deprecated) 2724* = 5: the minimum column dimension for blocking to be used; 2725* rectangular blocks must have dimension at least k by m, 2726* where k is given by ilaenv(2,...) and m by ilaenv(5,...) 2727* = 6: the crossover point for the svd (when reducing an m by n 2728* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds 2729* this value, a qr factorization is used first to reduce 2730* the matrix to a triangular form.) 2731* = 7: the number of processors 2732* = 8: the crossover point for the multishift qr method 2733* for nonsymmetric eigenvalue problems (deprecated) 2734* = 9: maximum size of the subproblems at the bottom of the 2735* computation tree in the divide-and-conquer algorithm 2736* (used by xgelsd and xgesdd) 2737* =10: ieee nan arithmetic can be trusted not to trap 2738* =11: infinity arithmetic can be trusted not to trap 2739* 12 <= ispec <= 16: 2740* xhseqr or one of its subroutines, 2741* see iparmq for detailed explanation 2742* 2743* name (input) character*(*) 2744* the name of the calling subroutine, in either upper case or 2745* lower case. 2746* 2747* opts (input) character*(*) 2748* the character options to the subroutine name, concatenated 2749* into a single character string. for example, uplo = 'u', 2750* trans = 't', and diag = 'n' for a triangular routine would 2751* be specified as opts = 'utn'. 2752* 2753* n1 (input) integer 2754* n2 (input) integer 2755* n3 (input) integer 2756* n4 (input) integer 2757* problem dimensions for the subroutine name; these may not all 2758* be required. 2759* 2760* further details 2761* =============== 2762* 2763* the following conventions have been used when calling ilaenv from the 2764* lapack routines: 2765* 1) opts is a concatenation of all of the character options to 2766* subroutine name, in the same order that they appear in the 2767* argument list for name, even if they are not used in determining 2768* the value of the parameter specified by ispec. 2769* 2) the problem dimensions n1, n2, n3, n4 are specified in the order 2770* that they appear in the argument list for name. n1 is used 2771* first, n2 second, and so on, and unused problem dimensions are 2772* passed a value of -1. 2773* 3) the parameter value returned by ilaenv is checked for validity in 2774* the calling subroutine. for example, ilaenv is used to retrieve 2775* the optimal blocksize for strtri as follows: 2776* 2777* nb = ilaenv( 1, 'strtri', uplo // diag, n, -1, -1, -1 ) 2778* if( nb.le.1 ) nb = max( 1, n ) 2779* 2780* ===================================================================== 2781* 2782* .. local scalars .. 2783 integer i, ic, iz, nb, nbmin, nx 2784 logical cname, sname 2785 character c1*1, c2*2, c4*2, c3*3, subnam*6 2786* .. 2787* .. intrinsic functions .. 2788 intrinsic char, ichar, int, min, real 2789* .. 2790* .. external functions .. 2791 integer ieeeck, iparmq 2792 external ieeeck, iparmq 2793* .. 2794* .. executable statements .. 2795* 2796 go to ( 10, 10, 10, 80, 90, 100, 110, 120, 2797 $ 130, 140, 150, 160, 160, 160, 160, 160 )ispec 2798* 2799* invalid value for ispec 2800* 2801 ilaenv = -1 2802 return 2803* 2804 10 continue 2805* 2806* convert name to upper case if the first character is lower case. 2807* 2808 ilaenv = 1 2809 subnam = name 2810 ic = ichar( subnam( 1: 1 ) ) 2811 iz = ichar( 'z' ) 2812 if( iz.eq.90 .or. iz.eq.122 ) then 2813* 2814* ascii character set 2815* 2816 if( ic.ge.97 .and. ic.le.122 ) then 2817 subnam( 1: 1 ) = char( ic-32 ) 2818 do 20 i = 2, 6 2819 ic = ichar( subnam( i: i ) ) 2820 if( ic.ge.97 .and. ic.le.122 ) 2821 $ subnam( i: i ) = char( ic-32 ) 2822 20 continue 2823 end if 2824* 2825 else if( iz.eq.233 .or. iz.eq.169 ) then 2826* 2827* ebcdic character set 2828* 2829 if( ( ic.ge.129 .and. ic.le.137 ) .or. 2830 $ ( ic.ge.145 .and. ic.le.153 ) .or. 2831 $ ( ic.ge.162 .and. ic.le.169 ) ) then 2832 subnam( 1: 1 ) = char( ic+64 ) 2833 do 30 i = 2, 6 2834 ic = ichar( subnam( i: i ) ) 2835 if( ( ic.ge.129 .and. ic.le.137 ) .or. 2836 $ ( ic.ge.145 .and. ic.le.153 ) .or. 2837 $ ( ic.ge.162 .and. ic.le.169 ) )subnam( i: 2838 $ i ) = char( ic+64 ) 2839 30 continue 2840 end if 2841* 2842 else if( iz.eq.218 .or. iz.eq.250 ) then 2843* 2844* prime machines: ascii+128 2845* 2846 if( ic.ge.225 .and. ic.le.250 ) then 2847 subnam( 1: 1 ) = char( ic-32 ) 2848 do 40 i = 2, 6 2849 ic = ichar( subnam( i: i ) ) 2850 if( ic.ge.225 .and. ic.le.250 ) 2851 $ subnam( i: i ) = char( ic-32 ) 2852 40 continue 2853 end if 2854 end if 2855* 2856 c1 = subnam( 1: 1 ) 2857 sname = c1.eq.'s' .or. c1.eq.'d' 2858 cname = c1.eq.'c' .or. c1.eq.'z' 2859 if( .not.( cname .or. sname ) ) 2860 $ return 2861 c2 = subnam( 2: 3 ) 2862 c3 = subnam( 4: 6 ) 2863 c4 = c3( 2: 3 ) 2864* 2865 go to ( 50, 60, 70 )ispec 2866* 2867 50 continue 2868* 2869* ispec = 1: block size 2870* 2871* in these examples, separate code is provided for setting nb for 2872* real and complex. we assume that nb will take the same value in 2873* single or double precision. 2874* 2875 nb = 1 2876* 2877 if( c2.eq.'ge' ) then 2878 if( c3.eq.'trf' ) then 2879 if( sname ) then 2880 nb = 64 2881 else 2882 nb = 64 2883 end if 2884 else if( c3.eq.'qrf' .or. c3.eq.'rqf' .or. c3.eq.'lqf' .or. 2885 $ c3.eq.'qlf' ) then 2886 if( sname ) then 2887 nb = 32 2888 else 2889 nb = 32 2890 end if 2891 else if( c3.eq.'hrd' ) then 2892 if( sname ) then 2893 nb = 32 2894 else 2895 nb = 32 2896 end if 2897 else if( c3.eq.'brd' ) then 2898 if( sname ) then 2899 nb = 32 2900 else 2901 nb = 32 2902 end if 2903 else if( c3.eq.'tri' ) then 2904 if( sname ) then 2905 nb = 64 2906 else 2907 nb = 64 2908 end if 2909 end if 2910 else if( c2.eq.'po' ) then 2911 if( c3.eq.'trf' ) then 2912 if( sname ) then 2913 nb = 64 2914 else 2915 nb = 64 2916 end if 2917 end if 2918 else if( c2.eq.'sy' ) then 2919 if( c3.eq.'trf' ) then 2920 if( sname ) then 2921 nb = 64 2922 else 2923 nb = 64 2924 end if 2925 else if( sname .and. c3.eq.'trd' ) then 2926 nb = 32 2927 else if( sname .and. c3.eq.'gst' ) then 2928 nb = 64 2929 end if 2930 else if( cname .and. c2.eq.'he' ) then 2931 if( c3.eq.'trf' ) then 2932 nb = 64 2933 else if( c3.eq.'trd' ) then 2934 nb = 32 2935 else if( c3.eq.'gst' ) then 2936 nb = 64 2937 end if 2938 else if( sname .and. c2.eq.'or' ) then 2939 if( c3( 1: 1 ).eq.'g' ) then 2940 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 2941 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 2942 $ then 2943 nb = 32 2944 end if 2945 else if( c3( 1: 1 ).eq.'m' ) then 2946 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 2947 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 2948 $ then 2949 nb = 32 2950 end if 2951 end if 2952 else if( cname .and. c2.eq.'un' ) then 2953 if( c3( 1: 1 ).eq.'g' ) then 2954 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 2955 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 2956 $ then 2957 nb = 32 2958 end if 2959 else if( c3( 1: 1 ).eq.'m' ) then 2960 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 2961 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 2962 $ then 2963 nb = 32 2964 end if 2965 end if 2966 else if( c2.eq.'gb' ) then 2967 if( c3.eq.'trf' ) then 2968 if( sname ) then 2969 if( n4.le.64 ) then 2970 nb = 1 2971 else 2972 nb = 32 2973 end if 2974 else 2975 if( n4.le.64 ) then 2976 nb = 1 2977 else 2978 nb = 32 2979 end if 2980 end if 2981 end if 2982 else if( c2.eq.'pb' ) then 2983 if( c3.eq.'trf' ) then 2984 if( sname ) then 2985 if( n2.le.64 ) then 2986 nb = 1 2987 else 2988 nb = 32 2989 end if 2990 else 2991 if( n2.le.64 ) then 2992 nb = 1 2993 else 2994 nb = 32 2995 end if 2996 end if 2997 end if 2998 else if( c2.eq.'tr' ) then 2999 if( c3.eq.'tri' ) then 3000 if( sname ) then 3001 nb = 64 3002 else 3003 nb = 64 3004 end if 3005 end if 3006 else if( c2.eq.'la' ) then 3007 if( c3.eq.'uum' ) then 3008 if( sname ) then 3009 nb = 64 3010 else 3011 nb = 64 3012 end if 3013 end if 3014 else if( sname .and. c2.eq.'st' ) then 3015 if( c3.eq.'ebz' ) then 3016 nb = 1 3017 end if 3018 end if 3019 ilaenv = nb 3020 return 3021* 3022 60 continue 3023* 3024* ispec = 2: minimum block size 3025* 3026 nbmin = 2 3027 if( c2.eq.'ge' ) then 3028 if( c3.eq.'qrf' .or. c3.eq.'rqf' .or. c3.eq.'lqf' .or. c3.eq. 3029 $ 'qlf' ) then 3030 if( sname ) then 3031 nbmin = 2 3032 else 3033 nbmin = 2 3034 end if 3035 else if( c3.eq.'hrd' ) then 3036 if( sname ) then 3037 nbmin = 2 3038 else 3039 nbmin = 2 3040 end if 3041 else if( c3.eq.'brd' ) then 3042 if( sname ) then 3043 nbmin = 2 3044 else 3045 nbmin = 2 3046 end if 3047 else if( c3.eq.'tri' ) then 3048 if( sname ) then 3049 nbmin = 2 3050 else 3051 nbmin = 2 3052 end if 3053 end if 3054 else if( c2.eq.'sy' ) then 3055 if( c3.eq.'trf' ) then 3056 if( sname ) then 3057 nbmin = 8 3058 else 3059 nbmin = 8 3060 end if 3061 else if( sname .and. c3.eq.'trd' ) then 3062 nbmin = 2 3063 end if 3064 else if( cname .and. c2.eq.'he' ) then 3065 if( c3.eq.'trd' ) then 3066 nbmin = 2 3067 end if 3068 else if( sname .and. c2.eq.'or' ) then 3069 if( c3( 1: 1 ).eq.'g' ) then 3070 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 3071 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 3072 $ then 3073 nbmin = 2 3074 end if 3075 else if( c3( 1: 1 ).eq.'m' ) then 3076 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 3077 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 3078 $ then 3079 nbmin = 2 3080 end if 3081 end if 3082 else if( cname .and. c2.eq.'un' ) then 3083 if( c3( 1: 1 ).eq.'g' ) then 3084 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 3085 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 3086 $ then 3087 nbmin = 2 3088 end if 3089 else if( c3( 1: 1 ).eq.'m' ) then 3090 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 3091 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 3092 $ then 3093 nbmin = 2 3094 end if 3095 end if 3096 end if 3097 ilaenv = nbmin 3098 return 3099* 3100 70 continue 3101* 3102* ispec = 3: crossover point 3103* 3104 nx = 0 3105 if( c2.eq.'ge' ) then 3106 if( c3.eq.'qrf' .or. c3.eq.'rqf' .or. c3.eq.'lqf' .or. c3.eq. 3107 $ 'qlf' ) then 3108 if( sname ) then 3109 nx = 128 3110 else 3111 nx = 128 3112 end if 3113 else if( c3.eq.'hrd' ) then 3114 if( sname ) then 3115 nx = 128 3116 else 3117 nx = 128 3118 end if 3119 else if( c3.eq.'brd' ) then 3120 if( sname ) then 3121 nx = 128 3122 else 3123 nx = 128 3124 end if 3125 end if 3126 else if( c2.eq.'sy' ) then 3127 if( sname .and. c3.eq.'trd' ) then 3128 nx = 32 3129 end if 3130 else if( cname .and. c2.eq.'he' ) then 3131 if( c3.eq.'trd' ) then 3132 nx = 32 3133 end if 3134 else if( sname .and. c2.eq.'or' ) then 3135 if( c3( 1: 1 ).eq.'g' ) then 3136 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 3137 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 3138 $ then 3139 nx = 128 3140 end if 3141 end if 3142 else if( cname .and. c2.eq.'un' ) then 3143 if( c3( 1: 1 ).eq.'g' ) then 3144 if( c4.eq.'qr' .or. c4.eq.'rq' .or. c4.eq.'lq' .or. c4.eq. 3145 $ 'ql' .or. c4.eq.'hr' .or. c4.eq.'tr' .or. c4.eq.'br' ) 3146 $ then 3147 nx = 128 3148 end if 3149 end if 3150 end if 3151 ilaenv = nx 3152 return 3153* 3154 80 continue 3155* 3156* ispec = 4: number of shifts (used by xhseqr) 3157* 3158 ilaenv = 6 3159 return 3160* 3161 90 continue 3162* 3163* ispec = 5: minimum column dimension (not used) 3164* 3165 ilaenv = 2 3166 return 3167* 3168 100 continue 3169* 3170* ispec = 6: crossover point for svd (used by xgelss and xgesvd) 3171* 3172 ilaenv = int( real( min( n1, n2 ) )*1.6e0 ) 3173 return 3174* 3175 110 continue 3176* 3177* ispec = 7: number of processors (not used) 3178* 3179 ilaenv = 1 3180 return 3181* 3182 120 continue 3183* 3184* ispec = 8: crossover point for multishift (used by xhseqr) 3185* 3186 ilaenv = 50 3187 return 3188* 3189 130 continue 3190* 3191* ispec = 9: maximum size of the subproblems at the bottom of the 3192* computation tree in the divide-and-conquer algorithm 3193* (used by xgelsd and xgesdd) 3194* 3195 ilaenv = 25 3196 return 3197* 3198 140 continue 3199* 3200* ispec = 10: ieee nan arithmetic can be trusted not to trap 3201* 3202* ilaenv = 0 3203 ilaenv = 1 3204 if( ilaenv.eq.1 ) then 3205 ilaenv = ieeeck( 1, 0.0, 1.0 ) 3206 end if 3207 return 3208* 3209 150 continue 3210* 3211* ispec = 11: infinity arithmetic can be trusted not to trap 3212* 3213* ilaenv = 0 3214 ilaenv = 1 3215 if( ilaenv.eq.1 ) then 3216 ilaenv = ieeeck( 0, 0.0, 1.0 ) 3217 end if 3218 return 3219* 3220 160 continue 3221* 3222* 12 <= ispec <= 16: xhseqr or one of its subroutines. 3223* 3224 ilaenv = iparmq( ispec, name, opts, n1, n2, n3, n4 ) 3225 return 3226* 3227* end of ilaenv 3228* 3229 end 3230ccc 3231 integer function ieeeck( ispec, zero, one ) 3232* 3233* -- lapack auxiliary routine (version 3.2.2) -- 3234* -- lapack is a software package provided by univ. of tennessee, -- 3235* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 3236* june 2010 3237* 3238* .. scalar arguments .. 3239 integer ispec 3240 real one, zero 3241* .. 3242* 3243* purpose 3244* ======= 3245* 3246* ieeeck is called from the ilaenv to verify that infinity and 3247* possibly nan arithmetic is safe (i.e. will not trap). 3248* 3249* arguments 3250* ========= 3251* 3252* ispec (input) integer 3253* specifies whether to test just for inifinity arithmetic 3254* or whether to test for infinity and nan arithmetic. 3255* = 0: verify infinity arithmetic only. 3256* = 1: verify infinity and nan arithmetic. 3257* 3258* zero (input) real 3259* must contain the value 0.0 3260* this is passed to prevent the compiler from optimizing 3261* away this code. 3262* 3263* one (input) real 3264* must contain the value 1.0 3265* this is passed to prevent the compiler from optimizing 3266* away this code. 3267* 3268* return value: integer 3269* = 0: arithmetic failed to produce the correct answers 3270* = 1: arithmetic produced the correct answers 3271* 3272* .. local scalars .. 3273 real nan1, nan2, nan3, nan4, nan5, nan6, neginf, 3274 $ negzro, newzro, posinf 3275* .. 3276* .. executable statements .. 3277 ieeeck = 1 3278* 3279 posinf = one / zero 3280 if( posinf.le.one ) then 3281 ieeeck = 0 3282 return 3283 end if 3284* 3285 neginf = -one / zero 3286 if( neginf.ge.zero ) then 3287 ieeeck = 0 3288 return 3289 end if 3290* 3291 negzro = one / ( neginf+one ) 3292 if( negzro.ne.zero ) then 3293 ieeeck = 0 3294 return 3295 end if 3296* 3297 neginf = one / negzro 3298 if( neginf.ge.zero ) then 3299 ieeeck = 0 3300 return 3301 end if 3302* 3303 newzro = negzro + zero 3304 if( newzro.ne.zero ) then 3305 ieeeck = 0 3306 return 3307 end if 3308* 3309 posinf = one / newzro 3310 if( posinf.le.one ) then 3311 ieeeck = 0 3312 return 3313 end if 3314* 3315 neginf = neginf*posinf 3316 if( neginf.ge.zero ) then 3317 ieeeck = 0 3318 return 3319 end if 3320* 3321 posinf = posinf*posinf 3322 if( posinf.le.one ) then 3323 ieeeck = 0 3324 return 3325 end if 3326* 3327* 3328* 3329* 3330* return if we were only asked to check infinity arithmetic 3331* 3332 if( ispec.eq.0 ) 3333 $ return 3334* 3335 nan1 = posinf + neginf 3336* 3337 nan2 = posinf / neginf 3338* 3339 nan3 = posinf / posinf 3340* 3341 nan4 = posinf*zero 3342* 3343 nan5 = neginf*negzro 3344* 3345 nan6 = nan5*zero 3346* 3347 if( nan1.eq.nan1 ) then 3348 ieeeck = 0 3349 return 3350 end if 3351* 3352 if( nan2.eq.nan2 ) then 3353 ieeeck = 0 3354 return 3355 end if 3356* 3357 if( nan3.eq.nan3 ) then 3358 ieeeck = 0 3359 return 3360 end if 3361* 3362 if( nan4.eq.nan4 ) then 3363 ieeeck = 0 3364 return 3365 end if 3366* 3367 if( nan5.eq.nan5 ) then 3368 ieeeck = 0 3369 return 3370 end if 3371* 3372 if( nan6.eq.nan6 ) then 3373 ieeeck = 0 3374 return 3375 end if 3376* 3377 return 3378 end 3379ccc 3380 integer function idamax(n,dx,incx) 3381* .. scalar arguments .. 3382 integer incx,n 3383* .. 3384* .. array arguments .. 3385 double precision dx(*) 3386* .. 3387* 3388* purpose 3389* ======= 3390* 3391* idamax finds the index of element having max. absolute value. 3392* 3393* further details 3394* =============== 3395* 3396* jack dongarra, linpack, 3/11/78. 3397* modified 3/93 to return if incx .le. 0. 3398* modified 12/3/93, array(1) declarations changed to array(*) 3399* 3400* ===================================================================== 3401* 3402* .. local scalars .. 3403 double precision dmax 3404 integer i,ix 3405* .. 3406* .. intrinsic functions .. 3407 intrinsic dabs 3408* .. 3409 idamax = 0 3410 if (n.lt.1 .or. incx.le.0) return 3411 idamax = 1 3412 if (n.eq.1) return 3413 if (incx.eq.1) go to 20 3414* 3415* code for increment not equal to 1 3416* 3417 ix = 1 3418 dmax = dabs(dx(1)) 3419 ix = ix + incx 3420 do 10 i = 2,n 3421 if (dabs(dx(ix)).le.dmax) go to 5 3422 idamax = i 3423 dmax = dabs(dx(ix)) 3424 5 ix = ix + incx 3425 10 continue 3426 return 3427* 3428* code for increment equal to 1 3429* 3430 20 dmax = dabs(dx(1)) 3431 do 30 i = 2,n 3432 if (dabs(dx(i)).le.dmax) go to 30 3433 idamax = i 3434 dmax = dabs(dx(i)) 3435 30 continue 3436 return 3437 end 3438ccc 3439 logical function lsame(ca,cb) 3440* 3441* -- lapack auxiliary routine (version 3.1) -- 3442* univ. of tennessee, univ. of california berkeley and nag ltd.. 3443* november 2006 3444* 3445* .. scalar arguments .. 3446 character ca,cb 3447* .. 3448* 3449* purpose 3450* ======= 3451* 3452* lsame returns .true. if ca is the same letter as cb regardless of 3453* case. 3454* 3455* arguments 3456* ========= 3457* 3458* ca (input) character*1 3459* 3460* cb (input) character*1 3461* ca and cb specify the single characters to be compared. 3462* 3463* ===================================================================== 3464* 3465* .. intrinsic functions .. 3466 intrinsic ichar 3467* .. 3468* .. local scalars .. 3469 integer inta,intb,zcode 3470* .. 3471* 3472* test if the characters are equal 3473* 3474 lsame = ca .eq. cb 3475 if (lsame) return 3476* 3477* now test for equivalence if both characters are alphabetic. 3478* 3479 zcode = ichar('z') 3480* 3481* use 'z' rather than 'a' so that ascii can be detected on prime 3482* machines, on which ichar returns a value with bit 8 set. 3483* ichar('a') on prime machines returns 193 which is the same as 3484* ichar('a') on an ebcdic machine. 3485* 3486 inta = ichar(ca) 3487 intb = ichar(cb) 3488* 3489 if (zcode.eq.90 .or. zcode.eq.122) then 3490* 3491* ascii is assumed - zcode is the ascii code of either lower or 3492* upper case 'z'. 3493* 3494 if (inta.ge.97 .and. inta.le.122) inta = inta - 32 3495 if (intb.ge.97 .and. intb.le.122) intb = intb - 32 3496* 3497 else if (zcode.eq.233 .or. zcode.eq.169) then 3498* 3499* ebcdic is assumed - zcode is the ebcdic code of either lower or 3500* upper case 'z'. 3501* 3502 if (inta.ge.129 .and. inta.le.137 .or. 3503 + inta.ge.145 .and. inta.le.153 .or. 3504 + inta.ge.162 .and. inta.le.169) inta = inta + 64 3505 if (intb.ge.129 .and. intb.le.137 .or. 3506 + intb.ge.145 .and. intb.le.153 .or. 3507 + intb.ge.162 .and. intb.le.169) intb = intb + 64 3508* 3509 else if (zcode.eq.218 .or. zcode.eq.250) then 3510* 3511* ascii is assumed, on prime machines - zcode is the ascii code 3512* plus 128 of either lower or upper case 'z'. 3513* 3514 if (inta.ge.225 .and. inta.le.250) inta = inta - 32 3515 if (intb.ge.225 .and. intb.le.250) intb = intb - 32 3516 end if 3517 lsame = inta .eq. intb 3518* 3519* return 3520* 3521* end of lsame 3522* 3523 end 3524ccc 3525 integer function iparmq( ispec, name, opts, n, ilo, ihi, lwork ) 3526* 3527* -- lapack auxiliary routine (version 3.2) -- 3528* -- lapack is a software package provided by univ. of tennessee, -- 3529* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 3530* november 2006 3531* 3532* .. scalar arguments .. 3533 integer ihi, ilo, ispec, lwork, n 3534 character name*( * ), opts*( * ) 3535* 3536* purpose 3537* ======= 3538* 3539* this program sets problem and machine dependent parameters 3540* useful for xhseqr and its subroutines. it is called whenever 3541* ilaenv is called with 12 <= ispec <= 16 3542* 3543* arguments 3544* ========= 3545* 3546* ispec (input) integer scalar 3547* ispec specifies which tunable parameter iparmq should 3548* return. 3549* 3550* ispec=12: (inmin) matrices of order nmin or less 3551* are sent directly to xlahqr, the implicit 3552* double shift qr algorithm. nmin must be 3553* at least 11. 3554* 3555* ispec=13: (inwin) size of the deflation window. 3556* this is best set greater than or equal to 3557* the number of simultaneous shifts ns. 3558* larger matrices benefit from larger deflation 3559* windows. 3560* 3561* ispec=14: (inibl) determines when to stop nibbling and 3562* invest in an (expensive) multi-shift qr sweep. 3563* if the aggressive early deflation subroutine 3564* finds ld converged eigenvalues from an order 3565* nw deflation window and ld.gt.(nw*nibble)/100, 3566* then the next qr sweep is skipped and early 3567* deflation is applied immediately to the 3568* remaining active diagonal block. setting 3569* iparmq(ispec=14) = 0 causes ttqre to skip a 3570* multi-shift qr sweep whenever early deflation 3571* finds a converged eigenvalue. setting 3572* iparmq(ispec=14) greater than or equal to 100 3573* prevents ttqre from skipping a multi-shift 3574* qr sweep. 3575* 3576* ispec=15: (nshfts) the number of simultaneous shifts in 3577* a multi-shift qr iteration. 3578* 3579* ispec=16: (iacc22) iparmq is set to 0, 1 or 2 with the 3580* following meanings. 3581* 0: during the multi-shift qr sweep, 3582* xlaqr5 does not accumulate reflections and 3583* does not use matrix-matrix multiply to 3584* update the far-from-diagonal matrix 3585* entries. 3586* 1: during the multi-shift qr sweep, 3587* xlaqr5 and/or xlaqraccumulates reflections and uses 3588* matrix-matrix multiply to update the 3589* far-from-diagonal matrix entries. 3590* 2: during the multi-shift qr sweep. 3591* xlaqr5 accumulates reflections and takes 3592* advantage of 2-by-2 block structure during 3593* matrix-matrix multiplies. 3594* (if xtrmm is slower than xgemm, then 3595* iparmq(ispec=16)=1 may be more efficient than 3596* iparmq(ispec=16)=2 despite the greater level of 3597* arithmetic work implied by the latter choice.) 3598* 3599* name (input) character string 3600* name of the calling subroutine 3601* 3602* opts (input) character string 3603* this is a concatenation of the string arguments to 3604* ttqre. 3605* 3606* n (input) integer scalar 3607* n is the order of the hessenberg matrix h. 3608* 3609* ilo (input) integer 3610* ihi (input) integer 3611* it is assumed that h is already upper triangular 3612* in rows and columns 1:ilo-1 and ihi+1:n. 3613* 3614* lwork (input) integer scalar 3615* the amount of workspace available. 3616* 3617* further details 3618* =============== 3619* 3620* little is known about how best to choose these parameters. 3621* it is possible to use different values of the parameters 3622* for each of chseqr, dhseqr, shseqr and zhseqr. 3623* 3624* it is probably best to choose different parameters for 3625* different matrices and different parameters at different 3626* times during the iteration, but this has not been 3627* implemented --- yet. 3628* 3629* 3630* the best choices of most of the parameters depend 3631* in an ill-understood way on the relative execution 3632* rate of xlaqr3 and xlaqr5 and on the nature of each 3633* particular eigenvalue problem. experiment may be the 3634* only practical way to determine which choices are most 3635* effective. 3636* 3637* following is a list of default values supplied by iparmq. 3638* these defaults may be adjusted in order to attain better 3639* performance in any particular computational environment. 3640* 3641* iparmq(ispec=12) the xlahqr vs xlaqr0 crossover point. 3642* default: 75. (must be at least 11.) 3643* 3644* iparmq(ispec=13) recommended deflation window size. 3645* this depends on ilo, ihi and ns, the 3646* number of simultaneous shifts returned 3647* by iparmq(ispec=15). the default for 3648* (ihi-ilo+1).le.500 is ns. the default 3649* for (ihi-ilo+1).gt.500 is 3*ns/2. 3650* 3651* iparmq(ispec=14) nibble crossover point. default: 14. 3652* 3653* iparmq(ispec=15) number of simultaneous shifts, ns. 3654* a multi-shift qr iteration. 3655* 3656* if ihi-ilo+1 is ... 3657* 3658* greater than ...but less ... the 3659* or equal to ... than default is 3660* 3661* 0 30 ns = 2+ 3662* 30 60 ns = 4+ 3663* 60 150 ns = 10 3664* 150 590 ns = ** 3665* 590 3000 ns = 64 3666* 3000 6000 ns = 128 3667* 6000 infinity ns = 256 3668* 3669* (+) by default matrices of this order are 3670* passed to the implicit double shift routine 3671* xlahqr. see iparmq(ispec=12) above. these 3672* values of ns are used only in case of a rare 3673* xlahqr failure. 3674* 3675* (**) the asterisks (**) indicate an ad-hoc 3676* function increasing from 10 to 64. 3677* 3678* iparmq(ispec=16) select structured matrix multiply. 3679* (see ispec=16 above for details.) 3680* default: 3. 3681* 3682* ================================================================ 3683* .. parameters .. 3684 integer inmin, inwin, inibl, ishfts, iacc22 3685 parameter ( inmin = 12, inwin = 13, inibl = 14, 3686 $ ishfts = 15, iacc22 = 16 ) 3687 integer nmin, k22min, kacmin, nibble, knwswp 3688 parameter ( nmin = 75, k22min = 14, kacmin = 14, 3689 $ nibble = 14, knwswp = 500 ) 3690 real two 3691 parameter ( two = 2.0 ) 3692* .. 3693* .. local scalars .. 3694 integer nh, ns 3695* .. 3696* .. intrinsic functions .. 3697 intrinsic log, max, mod, nint, real 3698* .. 3699* .. executable statements .. 3700 if( ( ispec.eq.ishfts ) .or. ( ispec.eq.inwin ) .or. 3701 $ ( ispec.eq.iacc22 ) ) then 3702* 3703* ==== set the number simultaneous shifts ==== 3704* 3705 nh = ihi - ilo + 1 3706 ns = 2 3707 if( nh.ge.30 ) 3708 $ ns = 4 3709 if( nh.ge.60 ) 3710 $ ns = 10 3711 if( nh.ge.150 ) 3712 $ ns = max( 10, nh / nint( log( real( nh ) ) / log( two ) ) ) 3713 if( nh.ge.590 ) 3714 $ ns = 64 3715 if( nh.ge.3000 ) 3716 $ ns = 128 3717 if( nh.ge.6000 ) 3718 $ ns = 256 3719 ns = max( 2, ns-mod( ns, 2 ) ) 3720 end if 3721* 3722 if( ispec.eq.inmin ) then 3723* 3724* 3725* ===== matrices of order smaller than nmin get sent 3726* . to xlahqr, the classic double shift algorithm. 3727* . this must be at least 11. ==== 3728* 3729 iparmq = nmin 3730* 3731 else if( ispec.eq.inibl ) then 3732* 3733* ==== inibl: skip a multi-shift qr iteration and 3734* . whenever aggressive early deflation finds 3735* . at least (nibble*(window size)/100) deflations. ==== 3736* 3737 iparmq = nibble 3738* 3739 else if( ispec.eq.ishfts ) then 3740* 3741* ==== nshfts: the number of simultaneous shifts ===== 3742* 3743 iparmq = ns 3744* 3745 else if( ispec.eq.inwin ) then 3746* 3747* ==== nw: deflation window size. ==== 3748* 3749 if( nh.le.knwswp ) then 3750 iparmq = ns 3751 else 3752 iparmq = 3*ns / 2 3753 end if 3754* 3755 else if( ispec.eq.iacc22 ) then 3756* 3757* ==== iacc22: whether to accumulate reflections 3758* . before updating the far-from-diagonal elements 3759* . and whether to use 2-by-2 block structure while 3760* . doing it. a small amount of work could be saved 3761* . by making this choice dependent also upon the 3762* . nh=ihi-ilo+1. 3763* 3764 iparmq = 0 3765 if( ns.ge.kacmin ) 3766 $ iparmq = 1 3767 if( ns.ge.k22min ) 3768 $ iparmq = 2 3769* 3770 else 3771* ===== invalid value of ispec ===== 3772 iparmq = -1 3773* 3774 end if 3775* 3776* ==== end of iparmq ==== 3777* 3778 end 3779ccc 3780 double precision function dlamch( cmach ) 3781* 3782* -- lapack auxiliary routine (version 3.3.0) -- 3783* -- lapack is a software package provided by univ. of tennessee, -- 3784* -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- 3785* based on lapack dlamch but with fortran 95 query functions 3786* see: http://www.cs.utk.edu/~luszczek/lapack/lamch.html 3787* and http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289 3788* july 2010 3789* 3790* .. scalar arguments .. 3791 character cmach 3792* .. 3793* 3794* purpose 3795* ======= 3796* 3797* dlamch determines double precision machine parameters. 3798* 3799* arguments 3800* ========= 3801* 3802* cmach (input) character*1 3803* specifies the value to be returned by dlamch: 3804* = 'e' or 'e', dlamch := eps 3805* = 's' or 's , dlamch := sfmin 3806* = 'b' or 'b', dlamch := base 3807* = 'p' or 'p', dlamch := eps*base 3808* = 'n' or 'n', dlamch := t 3809* = 'r' or 'r', dlamch := rnd 3810* = 'm' or 'm', dlamch := emin 3811* = 'u' or 'u', dlamch := rmin 3812* = 'l' or 'l', dlamch := emax 3813* = 'o' or 'o', dlamch := rmax 3814* 3815* where 3816* 3817* eps = relative machine precision 3818* sfmin = safe minimum, such that 1/sfmin does not overflow 3819* base = base of the machine 3820* prec = eps*base 3821* t = number of (base) digits in the mantissa 3822* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise 3823* emin = minimum exponent before (gradual) underflow 3824* rmin = underflow threshold - base**(emin-1) 3825* emax = largest exponent before overflow 3826* rmax = overflow threshold - (base**emax)*(1-eps) 3827* 3828* ===================================================================== 3829* 3830* .. parameters .. 3831 double precision one, zero 3832 parameter ( one = 1.0d+0, zero = 0.0d+0 ) 3833* .. 3834* .. local scalars .. 3835 double precision rnd, eps, sfmin, small, rmach 3836* .. 3837* .. external functions .. 3838 logical lsame 3839 external lsame 3840* .. 3841* .. intrinsic functions .. 3842 intrinsic digits, epsilon, huge, maxexponent, 3843 $ minexponent, radix, tiny 3844* .. 3845* .. executable statements .. 3846* 3847* 3848* assume rounding, not chopping. always. 3849* 3850 rnd = one 3851* 3852 if( one.eq.rnd ) then 3853 eps = epsilon(zero) * 0.5 3854 else 3855 eps = epsilon(zero) 3856 end if 3857* 3858 if( lsame( cmach, 'e' ) ) then 3859 rmach = eps 3860 else if( lsame( cmach, 's' ) ) then 3861 sfmin = tiny(zero) 3862 small = one / huge(zero) 3863 if( small.ge.sfmin ) then 3864* 3865* use small plus a bit, to avoid the possibility of rounding 3866* causing overflow when computing 1/sfmin. 3867* 3868 sfmin = small*( one+eps ) 3869 end if 3870 rmach = sfmin 3871 else if( lsame( cmach, 'b' ) ) then 3872 rmach = radix(zero) 3873 else if( lsame( cmach, 'p' ) ) then 3874 rmach = eps * radix(zero) 3875 else if( lsame( cmach, 'n' ) ) then 3876 rmach = digits(zero) 3877 else if( lsame( cmach, 'r' ) ) then 3878 rmach = rnd 3879 else if( lsame( cmach, 'm' ) ) then 3880 rmach = minexponent(zero) 3881 else if( lsame( cmach, 'u' ) ) then 3882 rmach = tiny(zero) 3883 else if( lsame( cmach, 'l' ) ) then 3884 rmach = maxexponent(zero) 3885 else if( lsame( cmach, 'o' ) ) then 3886 rmach = huge(zero) 3887 else 3888 rmach = zero 3889 end if 3890* 3891 dlamch = rmach 3892 return 3893* 3894* end of dlamch 3895* 3896 end 3897************************************************************************ 3898* 3899 double precision function dlamc3( a, b ) 3900* 3901* -- lapack auxiliary routine (version 3.3.0) -- 3902* univ. of tennessee, univ. of california berkeley and nag ltd.. 3903* november 2010 3904* 3905* .. scalar arguments .. 3906 double precision a, b 3907* .. 3908* 3909* purpose 3910* ======= 3911* 3912* dlamc3 is intended to force a and b to be stored prior to doing 3913* the addition of a and b , for use in situations where optimizers 3914* might hold one of these in a register. 3915* 3916* arguments 3917* ========= 3918* 3919* a (input) double precision 3920* b (input) double precision 3921* the values a and b. 3922* 3923* ===================================================================== 3924* 3925* .. executable statements .. 3926* 3927 dlamc3 = a + b 3928* 3929 return 3930* 3931* end of dlamc3 3932* 3933 end 3934* 3935************************************************************************ 3936c $Id: rism1d.F 21176 2011-10-10 06:35:49Z d3y133 $ 3937