1 subroutine mpolefit(idip,nesp,esp,connl,dx,dy,dz,iz,dmachg, 2 & keywrd,ichadd) 3 4c Routine for fitting of atomic monopoles, dipoles and quadrupoles to 5c the electrostatic potential. Multipoles are obtained in spherical 6c tensor form, in coordinate system of the molecule, 7c and are given in atomic units. 8c 9c Q_00 10c Q_10 Q_1c Q_1s 11c Q_20 Q_21c Q_21s Q_22c Q_22s 12c 13c Fitting strategy is based on: 14c Hinsen and Roux, J. Comp. Chem. 1997, 18, 368 15c using Singular Value Decomposition 16c constraints are imposed by elimination 17c 18c Wijnand Mooij, may 1998 19 20c 21c process keywords ATPOL MPTOL MPEQUIV 22c 23c ATPOL sets the use of charge, dipole and quadrupole on each atom 24c Default is M+D+Q, except for hydrogen M+D 25c 26c ATPOL=(atnum1/icode,atnum2-atnum10/icode,...) 27c 28c icode = +1 (if monopole) +2 (if dipole) +4 (if quadrupole) 29c 30c SV tolerance (default 1D-5) 31c 32c MPTOL=value 33c 34c Equivalences can be set for groups of atoms. 35c Note that equivalences are set in variable numbers, not in atom number; 36c these will differ when dipoles, quadrupoles are used. 37c examp: if atom 1=M+D+Q, then charge on atom 2 has variable number 10 38c 39c MPEQUIV=(var1,var2,var3/var4,var5/...) 40c 41 42 implicit double precision (a-h,o-z) 43 logical dmachg 44 dimension esp(*), connl(3,*) 45 parameter (numatm=2000) 46 parameter (mxel=100) 47 parameter (mesp = 30000) 48 parameter (np = 200) 49 parameter (tolc=1d-10) 50 parameter (tol=1d-5) 51 parameter (lfmax=2) 52 parameter (lfmaxf=(lfmax+1)*(lfmax+1)) 53 character*2 elemnt 54 common /elem/elemnt(mxel) 55 common /coord / xyz(3,numatm) 56 common /moldat/ natoms, norbs, nelecs,nat(numatm) 57 common /rdwr/ iun1,iun2,iun3,iun4,iun5 58 character*137 line 59 common /curlin/ line 60 character*(*) keywrd 61 character*137 keyhlp,tstr,strt 62 logical keyr 63 dimension u(mesp,np),v(np,np),w(np) 64 dimension b(np,np),c(np),p(np,np),binvc(np),esp1(mesp) 65 dimension scratch(np) 66 dimension pvec(3),temp(5) 67 dimension ict(50) 68 logical monopole,dipole,quadrupole 69 common /mpfit/ q(lfmaxf,numatm), 70 & monopole(numatm),dipole(numatm),quadrupole(numatm) 71 common /charge/ dipo(3),ihasq,ihsdp,iqon,idipon 72 73 write(iun3,*) ' ---------------- ' 74 write(iun3,*) ' MPOLEFIT' 75 write(iun3,*) ' ---------------- ' 76 77 do i=1,numatm 78 monopole(i) = .true. 79 dipole(i) = .false. 80 quadrupole(i) = .false. 81 end do 82 83 do i=natoms+1,numatm 84 nat(i) = 99 85 end do 86 87 elemnt(99) = ' X' 88 89 call setmp(keywrd) 90 91c Check for user-input SV tolerance 92 93 if (.not.keyr(keywrd,'MPTOL',toll)) then 94 toll = tol 95 endif 96 write(iun3,*) 'MPTOL =',toll 97 98 rt3 = dsqrt(3.0d0) 99c 100c conversion factor for debye to atomic units 101c 102 cf = 5.2917715d-11*1.601917d-19 / 3.33564d-30 103 toang = 0.52917706d0 104 105c calculate total charge on the molecule 106 107 iz = 0 108 do i=1,natoms 109 iz = iz + nat(i) 110 end do 111 iz = iz - nelecs 112 if (iz.ne.0) write(iun3,*) 'Charge of molecule = ',iz 113 114 if (ichadd.ne.0) then 115 116 call chadd(natoms) 117 118 do j=1,natoms 119 write(iun3,'(i4,3f10.5,i4,a)') j,(xyz(m,j)*toang,m=1,3), 120 & nat(j),elemnt(nat(j)) 121 end do 122 123 endif 124 125 126c Set the number of variables 127 128 nvar = 0 129 do i=1,natoms 130 if (monopole(i)) nvar = nvar + 1 131 if (dipole(i)) nvar = nvar + 3 132 if (quadrupole(i)) nvar = nvar + 5 133 end do 134 135 if (nvar.gt.np) then 136 write(iun3,*) 'Too many parameters in mpolefit' 137 stop 138 endif 139 if (nesp.gt.mesp) then 140 write(iun3,*) 'Too many ESP points in mpolefit' 141 stop 142 endif 143 144c Set up the constraint relation b*q=c 145 146 ivar = 0 147 do j=1,natoms 148 149 if (monopole(j)) then 150 ivar = ivar + 1 151 b(1,ivar) = 1.d0 152 if (idip.eq.1) then 153 b(2,ivar) = xyz(1,j) 154 b(3,ivar) = xyz(2,j) 155 b(4,ivar) = xyz(3,j) 156 endif 157 endif 158 159 if (dipole(j)) then 160 do k=1,3 161 b(1,ivar+k) = 0.d0 162 end do 163 if (idip.eq.1) then 164 165c X-component of dipole 166 167 b(2,ivar+1) = 0.d0 168 b(2,ivar+2) = 1.d0 169 b(2,ivar+3) = 0.d0 170 171c Y-component of dipole 172 173 b(3,ivar+1) = 0.d0 174 b(3,ivar+2) = 0.d0 175 b(3,ivar+3) = 1.d0 176 177c Z-component of dipole 178 179 b(4,ivar+1) = 1.d0 180 b(4,ivar+2) = 0.d0 181 b(4,ivar+3) = 0.d0 182 endif 183 ivar = ivar + 3 184 endif 185 186 if (quadrupole(j)) then 187 do k=1,5 188 b(1,ivar+k) = 0.d0 189 end do 190 if (idip.eq.1) then 191 do k=1,5 192 do i=2,4 193 b(i,ivar+k) = 0.d0 194 end do 195 end do 196 endif 197 ivar = ivar + 5 198 endif 199 200 end do 201 202 c(1) = dfloat(iz) 203 ncon = 1 204 205 if (idip.eq.1) then 206 c(2) = dx/cf 207 c(3) = dy/cf 208 c(4) = dz/cf 209 ncon = 4 210 endif 211 212c Augment the constraint matrix to square 213 214 do i=ncon+1,nvar 215 do j=1,nvar 216 b(i,j) = 0.d0 217 end do 218 c(i) = 0.d0 219 end do 220 221 222c Set equivalences in constraint matrix 223 224 i = index(keywrd,'MPEQUIV') 225 if (i.ne.0) then 226 i1 = index(keywrd(i+1:),'(') 227 if (i1.ne.0) then 228 i2 = index(keywrd(i+1:),')') 229 if (i2.ne.0) then 230 keyhlp = keywrd(i+i1+1:i+i2-1) 231 l = i2 - i1 - 1 232 call spatrm(keyhlp,l) 233 call setlin(keyhlp,47) 234 keyhlp = line(1:l) 235 do while (nxtwrd(tstr,nstr,itype,rtype).eq.1) 236 keyhlp = line(1:l) 237 call setlin(tstr,44) 238 n = 0 239 do while (nxtwrd(strt,nstr,itype,rtype).eq.2) 240 n = n + 1 241 ict(n) = itype 242 end do 243 do k=2,n 244 ncon = ncon + 1 245 b(ncon,ict(1)) = 1.d0 246 b(ncon,ict(k)) = -1.d0 247c print*,'b(',ncon,',',ict(1),') = 1.0d0', 248c & ' b(',ncon,',',ict(k),') = -1.0d0' 249 end do 250 call setlin(keyhlp,0) 251 end do 252 endif 253 endif 254 endif 255 256 13 format(6f9.4) 257 258 write(iun3,*) 'Number of variables : ',nvar 259 write(iun3,*) 'Number of constraints : ',ncon 260 261 if (ncon.gt.nvar) then 262 write(iun3,*) 'More constraints than variables in mpolefit' 263 stop 264 endif 265 266 267c Copy b to u. B is needed later on. 268 269 do i=1,nvar 270 do j=1,nvar 271 u(j,i) = b(j,i) 272 end do 273 end do 274 275 call svd(mesp,nvar,nvar,u,w,.true.,u,.true.,v,ierr,scratch,np) 276 277 wmax = 0.d0 278 do j=1,nvar 279 if (w(j).gt.wmax) wmax = w(j) 280 end do 281 282 283 thresh = tolc*wmax 284 if (thresh.lt.1d-10) thresh = 1d-10 285 286 nsvd = 0 287 do j=1,nvar 288 if (w(j).lt.thresh) then 289 nsvd = nsvd + 1 290 w(j) = 0.d0 291 endif 292 end do 293 write(iun3,*) 'Rank of constraint matrix : ',nvar-nsvd 294 295c Calculate 1/w(i)*Utranspose, in p 296 297 do i=1,nvar 298 if (w(i).ne.0.d0) then 299 do j=1,nvar 300 p(i,j) = u(j,i) / w(i) 301 end do 302 else 303 do j=1,nvar 304 p(i,j) = 0.d0 305 end do 306 endif 307 end do 308 309c Inverse of B, in u 310 311 do i=1,nvar 312 do j=1,nvar 313 u(i,j) = 0.d0 314 do k=1,nvar 315 u(i,j) = u(i,j) + v(i,k)*p(k,j) 316 end do 317 end do 318 end do 319 320c Calculate p=1-binv*b 321 322 do i=1,nvar 323 do j=1,nvar 324 p(i,j) = 0.d0 325 do k=1,nvar 326 p(i,j) = p(i,j) - u(i,k)*b(k,j) 327 end do 328 end do 329 p(i,i) = p(i,i) + 1.d0 330 end do 331 332c Calculate binv*c 333 334 do i=1,nvar 335 binvc(i) = 0.0d0 336 do j=1,nvar 337 binvc(i) = binvc(i) + u(i,j)*c(j) 338 end do 339 end do 340 341 342 do i=1,nesp 343 do j=1,nvar 344 u(i,j) = 0.d0 345 end do 346 esp1(i) = esp(i) 347 end do 348 349c Set up design matrix A`=AP en p`=p-A*Binv*C 350 351 do i=1,nesp 352 ivar = 0 353 do j=1,natoms 354 355 if (monopole(j)) then 356 ivar = ivar + 1 357 pvec(1) = connl(1,i) - xyz(1,j) 358 pvec(2) = connl(2,i) - xyz(2,j) 359 pvec(3) = connl(3,i) - xyz(3,j) 360 qp1 = pvec(1)*pvec(1) 361 qp2 = pvec(2)*pvec(2) 362 qp3 = pvec(3)*pvec(3) 363 r2 = qp1 + qp2 + qp3 364 r = dsqrt(r2) 365 temp1 = 1.d0 / r 366 do k=1,nvar 367 u(i,k) = u(i,k) + temp1*p(ivar,k) 368 end do 369 esp1(i) = esp1(i) - temp1*binvc(ivar) 370 endif 371 372 if (dipole(j)) then 373 pvec(1) = connl(1,i) - xyz(1,j) 374 pvec(2) = connl(2,i) - xyz(2,j) 375 pvec(3) = connl(3,i) - xyz(3,j) 376 qp1 = pvec(1)*pvec(1) 377 qp2 = pvec(2)*pvec(2) 378 qp3 = pvec(3)*pvec(3) 379 r2 = qp1 + qp2 + qp3 380 r = dsqrt(r2) 381 r3= r2*r 382 temp(1) = pvec(3) / r3 383 temp(2) = pvec(1) / r3 384 temp(3) = pvec(2) / r3 385 do kk=1,3 386 do k=1,nvar 387 u(i,k) = u(i,k) + temp(kk)*p(ivar+kk,k) 388 end do 389 esp1(i) = esp1(i) - temp(kk)*binvc(ivar+kk) 390 end do 391 ivar = ivar + 3 392 endif 393 394 if (quadrupole(j)) then 395 pvec(1) = connl(1,i) - xyz(1,j) 396 pvec(2) = connl(2,i) - xyz(2,j) 397 pvec(3) = connl(3,i) - xyz(3,j) 398 qp1 = pvec(1)*pvec(1) 399 qp2 = pvec(2)*pvec(2) 400 qp3 = pvec(3)*pvec(3) 401 r2 = qp1 + qp2 + qp3 402 r = dsqrt(r2) 403 r5 = r2*r2*r 404 temp(1) = 0.5d0 * (3.0d0*qp3 - r2)/ r5 405 temp(2) = rt3*pvec(1)*pvec(3)/r5 406 temp(3) = rt3*pvec(2)*pvec(3)/r5 407 temp(4) = 0.5d0*rt3*(qp1 - qp2)/r5 408 temp(5) = rt3*pvec(1)*pvec(2)/r5 409 do kk=1,5 410 do k=1,nvar 411 u(i,k) = u(i,k) + temp(kk)*p(ivar+kk,k) 412 end do 413 esp1(i) = esp1(i) - temp(kk)*binvc(ivar+kk) 414 end do 415 ivar = ivar + 5 416 endif 417 end do 418 end do 419 420 call svd(mesp,nesp,nvar,u,w,.true.,u,.true.,v,ierr,scratch,np) 421 422 wmax = 0.d0 423 do j=1,nvar 424 if (w(j).gt.wmax) wmax = w(j) 425 end do 426 thresh = toll*wmax 427 if (thresh.lt.1d-10) thresh = 1d-10 428 429 nsvd = 0 430 do j=1,nvar 431 if (w(j).lt.thresh) then 432 w(j) = 0.d0 433 nsvd = nsvd + 1 434 endif 435 end do 436 437 write(iun3,*) 'Rank of design matrix AP : ',nvar-nsvd 438 439c Solve for p` 440 441 do i=1,nvar 442 s = 0.d0 443 if (w(i).ne.0.d0) then 444 do j=1,nesp 445 s = s + u(j,i)*esp1(j) 446 end do 447 s = s / w(i) 448 endif 449 scratch(i) = s 450 end do 451 do i=1,nvar 452 s = 0.d0 453 do j=1,nvar 454 s = s + v(i,j)*scratch(j) 455 end do 456 c(i) = s 457 end do 458 459 460c Solution obeying constraints from q=BinvC + p*c 461c stored in w 462 463 do i=1,nvar 464 w(i) = binvc(i) 465 do j=1,nvar 466 w(i) = w(i) + p(i,j)*c(j) 467 end do 468 end do 469 470c Store fitted multipoles 471 472 ivar = 0 473 do i=1, natoms 474 475 if (monopole(i)) then 476 ivar = ivar + 1 477 q(1,i) = w(ivar) 478 else 479 q(1,i) = 0.d0 480 endif 481 482 if (dipole(i)) then 483 do j=1,3 484 ivar = ivar + 1 485 q(1+j,i) = w(ivar) 486 end do 487 else 488 do j=1,3 489 q(1+j,i) = 0.d0 490 end do 491 endif 492 493 if (quadrupole(i)) then 494 do j=1,5 495 ivar = ivar + 1 496 q(4+j,i) = w(ivar) 497 end do 498 else 499 do j=1,5 500 q(4+j,i) = 0.d0 501 end do 502 endif 503 504 end do 505 506c 507c calculate root mean square fits and relative root mean square fits 508c 509 rms = 0.d0 510 rrms = 0.d0 511 do i=1,nesp 512 espc = 0.d0 513 call calc2(connl(1,i),connl(2,i),connl(3,i),espc) 514 rms = rms + (espc-esp(i))**2 515 rrms = rrms + esp(i)**2 516 end do 517 518 rms = dsqrt(rms/nesp) 519 rrms = rms / dsqrt(rrms/nesp) 520 rms = rms*627.51d0 521 522 write(iun3,'(5x,''ATOM NO. TYPE'')') 523 write(iun3,*)' ' 524 525 qtot = 0.d0 526 do i=1,natoms 527 qtot = qtot + q(1,i) 528 write(iun3,'(7x,i2,9x,a2,1x,f9.5)')i,elemnt(nat(i)),q(1,i) 529 write(iun3,'(21x,3f9.5)') (q(j,i),j=2,4) 530 write(iun3,'(21x,5f9.5)') (q(j,i),j=5,9) 531 end do 532 533 write(iun3,*)' ' 534 write(iun3,*) 'THE TOTAL CHARGE IS: ',qtot 535 write(iun3,*)' ' 536 537 write(iun3,*) 'THE NUMBER OF POINTS IS: ',nesp 538 write(iun3,*) 'THE RMS DEVIATION IS: ',rms 539 write(iun3,*) 'THE RRMS DEVIATION IS: ',rrms 540 541 write(iun3,*)' ' 542 write(iun3,*) 543 & 'DIPOLE MOMENT EVALUATED FROM Monopoles and Dipoles (debye)' 544 write(iun3,'(12x,'' X Y Z TOTAL'')') 545 write(iun3,*)' ' 546 547 dipx = 0 548 dipy = 0 549 dipz = 0 550 551 do i=1,natoms 552 dipx = dipx + xyz(1,i)*q(1,i) 553 dipy = dipy + xyz(2,i)*q(1,i) 554 dipz = dipz + xyz(3,i)*q(1,i) 555 end do 556 557 dip = dsqrt(dipx**2+dipy**2+dipz**2) 558 559 write(iun3,'(8x,4f9.4,4x,"Charge")') 560 & dipx*cf,dipy*cf,dipz*cf,dip*cf 561 dipx = 0 562 dipy = 0 563 dipz = 0 564 565 do i=1,natoms 566 dipx = dipx + q(3,i) 567 dipy = dipy + q(4,i) 568 dipz = dipz + q(2,i) 569 end do 570 571 dip = dsqrt(dipx**2+dipy**2+dipz**2) 572 write(iun3,'(8x,4f9.4,4x,"Dipole")') 573 & dipx*cf,dipy*cf,dipz*cf,dip*cf 574 575 dipx = 0 576 dipy = 0 577 dipz = 0 578 579 do i=1,natoms 580 dipx = dipx + xyz(1,i)*q(1,i) 581 dipy = dipy + xyz(2,i)*q(1,i) 582 dipz = dipz + xyz(3,i)*q(1,i) 583 dipx = dipx + q(3,i) 584 dipy = dipy + q(4,i) 585 dipz = dipz + q(2,i) 586 end do 587 588 dip = dsqrt(dipx**2+dipy**2+dipz**2) 589 write(iun3,'(8x,4f9.4,4x,"Total/")') 590 & dipx*cf,dipy*cf,dipz*cf,dip*cf 591 592 dipo(1) = dipx 593 dipo(2) = dipy 594 dipo(3) = dipz 595 ihsdp = 3 596 597 call cpypol 598 599 open(unit=46,form='formatted',file='esp.xyz', 600 & status='unknown',err=200) 601 602 write(46,'(i5)') natoms 603 write(46,'(a)') 'Molden MPOLEFIT fitted charges' 604 do i=1,natoms 605 write(46,'(a2,1x,3(f9.6,1x),f9.6)') 606 & elemnt(nat(i)),(xyz(j,i)*0.52917706d0,j=1,3),q(1,i) 607 end do 608 close(46) 609 610 return 611 612 200 write(iun3,*) 'Couldnt write XYZ file esp.xyz' 613 end 614 615 616 subroutine setmp(keywrd) 617 implicit double precision (a-h,o-z) 618 parameter (numatm=2000) 619 parameter (lfmax=2) 620 parameter (lfmaxf=(lfmax+1)*(lfmax+1)) 621 logical monopole,dipole,quadrupole 622 common /mpfit/ q(lfmaxf,numatm), 623 & monopole(numatm),dipole(numatm),quadrupole(numatm) 624 common /moldat/ natoms, norbs, nelecs,nat(numatm) 625 character*(*) keywrd 626 dimension poles(numatm) 627 628 indx = index(keywrd,'ATPOL') 629 630 if (indx.ne.0) then 631 do i=1,natoms 632 monopole(i) = .false. 633 dipole(i) = .false. 634 quadrupole(i) = .false. 635 poles(i) = 7 636 if (nat(i).eq.1) poles(i) = 3 637 end do 638 call occin(indx,poles,numatm) 639 do i=1,natoms 640 j = int(poles(i)) 641 if (j.ge.4) then 642 quadrupole(i) = .true. 643 j = j - 4 644 endif 645 if (j.ge.2) then 646 dipole(i) = .true. 647 j = j - 2 648 endif 649 if (j.eq.1) monopole(i) = .true. 650 end do 651 else 652 do i=1,natoms 653 monopole(i) = .true. 654 dipole(i) = .true. 655 quadrupole(i) = .true. 656 if (nat(i).eq.1) then 657 quadrupole(i)=.false. 658 endif 659 end do 660 endif 661 662 return 663 end 664 665 subroutine calc2(x,y,z,pot) 666 implicit double precision (a-h,o-z), integer ( i-n) 667 parameter (lfmax=2) 668 parameter (lfmaxf=(lfmax+1)*(lfmax+1)) 669 parameter (numatm=2000) 670 logical monopole,dipole,quadrupole 671 common /mpfit/ q(lfmaxf,numatm), 672 & monopole(numatm),dipole(numatm),quadrupole(numatm) 673 common /coord/ xyz(3,numatm) 674 common /moldat/ natoms, norbs, nelecs,nat(numatm) 675 676 dimension pvec(3) 677c in this procedure the calculation order is inversed 678c we use 679c 680c a b c d 681c ---- + ---- + ---- + ---- = ((((d /x**2+c) /x**2)+b) /x**2) +a) /x 682c x x**3 x**5 x**7 683c 684 rt3=dsqrt(3.0d0) 685 686 pot=0.0d0 687 do 10 i=1,natoms 688 pvec(1)=x-xyz(1,i) 689 pvec(2)=y-xyz(2,i) 690 pvec(3)=z-xyz(3,i) 691 qp1=pvec(1)*pvec(1) 692 qp2=pvec(2)*pvec(2) 693 qp3=pvec(3)*pvec(3) 694 r2= qp1+qp2+qp3 695 r=dsqrt(r2) 696 pot = pot 697 & +((( 698 & (q(5,i)*0.5d0*(3.0d0*qp3-r2)+ 699 & q(6,i)*rt3*pvec(1)*pvec(3)+ 700 & q(7,i)*rt3*pvec(2)*pvec(3)+ 701 & q(8,i)*0.5d0*rt3*(qp1-qp2)+ 702 & q(9,i)*rt3*pvec(1)*pvec(2)))/r2 703 & +(q(2,i)*pvec(3)+q(3,i)*pvec(1)+q(4,i)*pvec(2)))/r2 704 & +q(1,i))/r 70510 continue 706 707 return 708 end 709 710 subroutine cpypol 711 implicit double precision (a-h,o-z), integer ( i-n) 712 parameter (numatm=2000) 713 parameter (lfmax=2) 714 parameter (lfmaxf=(lfmax+1)*(lfmax+1)) 715 parameter (mxsite=300) 716 parameter (lmax=4) 717 parameter (lmaxf=(lmax+1)*(lmax+1)) 718 logical monopole,dipole,quadrupole 719 common /mpfit/ q(lfmaxf,numatm), 720 & monopole(numatm),dipole(numatm),quadrupole(numatm) 721 character*8 ctag 722 common /multip/ qmom(lmaxf,mxsite),car(3,mxsite),ctag(mxsite), 723 & nsites 724 common /moldat/ natoms, norbs, nelecs,nat(numatm) 725 common /coord / xyz(3,numatm) 726 727 do i=1,natoms 728 do j=1,lmaxf 729 qmom(j,i) = 0.0d0 730 end do 731 do j=1,lfmaxf 732 qmom(j,i) = q(j,i) 733 end do 734 do j=1,3 735 car(j,i) = xyz(j,i) 736 end do 737 end do 738 739 nsites = natoms 740 741 return 742 end 743 744 745 SUBROUTINE SVD (NM, M, N, A, W, MATU, U, MATV, V, IERR, RV1, nn) 746C 747C***PURPOSE Perform the singular value decomposition of a rectangular 748C matrix. 749C***LIBRARY SLATEC 750C 751C This subroutine is a translation of the ALGOL procedure SVD, 752C NUM. MATH. 14, 403-420(1970) by Golub and Reinsch. 753C HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971). 754C 755C This subroutine determines the singular value decomposition 756C T 757C A=USV of a REAL M by N rectangular matrix. Householder 758C bidiagonalization and a variant of the QR algorithm are used. 759C 760C On Input 761C 762C NM must be set to the row dimension of the two-dimensional 763C array parameters, A, U, V, as declared in the calling 764C program dimension statement. NM is an INTEGER variable. 765C Note that NM must be at least as large as the maximum 766C of M and N. 767C 768C M is the number of rows of A and U. 769C 770C N is the number of columns of A and U and the order of V. 771C 772C A contains the rectangular input matrix to be decomposed. A is 773C a two-dimensional REAL array, dimensioned A(NM,N). 774C 775C MATU should be set to .TRUE. if the U matrix in the 776C decomposition is desired, and to .FALSE. otherwise. 777C MATU is a LOGICAL variable. 778C 779C MATV should be set to .TRUE. if the V matrix in the 780C decomposition is desired, and to .FALSE. otherwise. 781C MATV is a LOGICAL variable. 782C 783C On Output 784C 785C A is unaltered (unless overwritten by U or V). 786C 787C W contains the N (non-negative) singular values of A (the 788C diagonal elements of S). They are unordered. If an 789C error exit is made, the singular values should be correct 790C for indices IERR+1, IERR+2, ..., N. W is a one-dimensional 791C REAL array, dimensioned W(N). 792C 793C U contains the matrix U (orthogonal column vectors) of the 794C decomposition if MATU has been set to .TRUE. Otherwise, 795C U is used as a temporary array. U may coincide with A. 796C If an error exit is made, the columns of U corresponding 797C to indices of correct singular values should be correct. 798C U is a two-dimensional REAL array, dimensioned U(NM,N). 799C 800C V contains the matrix V (orthogonal) of the decomposition if 801C MATV has been set to .TRUE. Otherwise, V is not referenced. 802C V may also coincide with A if U does not. If an error 803C exit is made, the columns of V corresponding to indices of 804C correct singular values should be correct. V is a two- 805C dimensional REAL array, dimensioned V(NM,N). 806c+++ Changed: V can not coincide with A, as for memory reasons v is limited 807c to V(nn,nn). 808 809C 810C IERR is an INTEGER flag set to 811C Zero for normal return, 812C K if the K-th singular value has not been 813C determined after 30 iterations. 814C 815C RV1 is a one-dimensional REAL array used for temporary storage, 816C dimensioned RV1(N). 817c+++Changed: NN must be set to the row dimension of the two-dimensional 818C array parameters V, as declared in the calling 819C program dimension statement. NN is an INTEGER variable. 820C Note that NN must be at least as large as the maximum of N. 821 822c 823c 824 825C 826 827C 828C Questions and comments should be directed to B. S. Garbow, 829C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 830C ------------------------------------------------------------------ 831C 832 INTEGER I,J,K,L,M,N,II,I1,KK,K1,LL,L1,MN,NM,ITS,IERR,nn 833 double precision A(NM,*),W(*),U(NM,*),V(nn,*),RV1(*) 834 double precision C,F,G,H,S,X,Y,Z,SCALE,S1 835 LOGICAL MATU,MATV 836C 837 IERR = 0 838 839C 840 DO 100 I = 1, M 841C 842 DO 100 J = 1, N 843 U(I,J) = A(I,J) 844 100 CONTINUE 845C .......... HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM .......... 846 G = 0.0d0 847 SCALE = 0.0d0 848 S1 = 0.0d0 849C 850 DO 300 I = 1, N 851 L = I + 1 852 RV1(I) = SCALE * G 853 G = 0.0d0 854 S = 0.0d0 855 SCALE = 0.0d0 856 IF (I .GT. M) GO TO 210 857C 858 DO 120 K = I, M 859 120 SCALE = SCALE + ABS(U(K,I)) 860C 861 IF (SCALE .EQ. 0.0d0) GO TO 210 862C 863 DO 130 K = I, M 864 U(K,I) = U(K,I) / SCALE 865 S = S + U(K,I)**2 866 130 CONTINUE 867C 868 F = U(I,I) 869 G = -SIGN(SQRT(S),F) 870 H = F * G - S 871 U(I,I) = F - G 872 IF (I .EQ. N) GO TO 190 873C 874 DO 150 J = L, N 875 S = 0.0d0 876C 877 DO 140 K = I, M 878 140 S = S + U(K,I) * U(K,J) 879C 880 F = S / H 881C 882 DO 150 K = I, M 883 U(K,J) = U(K,J) + F * U(K,I) 884 150 CONTINUE 885C 886 190 DO 200 K = I, M 887 200 U(K,I) = SCALE * U(K,I) 888C 889 210 W(I) = SCALE * G 890 G = 0.0d0 891 S = 0.0d0 892 SCALE = 0.0d0 893 IF (I .GT. M .OR. I .EQ. N) GO TO 290 894C 895 DO 220 K = L, N 896 220 SCALE = SCALE + ABS(U(I,K)) 897C 898 IF (SCALE .EQ. 0.0d0) GO TO 290 899C 900 DO 230 K = L, N 901 U(I,K) = U(I,K) / SCALE 902 S = S + U(I,K)**2 903 230 CONTINUE 904C 905 F = U(I,L) 906 G = -SIGN(SQRT(S),F) 907 H = F * G - S 908 U(I,L) = F - G 909C 910 DO 240 K = L, N 911 240 RV1(K) = U(I,K) / H 912C 913 IF (I .EQ. M) GO TO 270 914C 915 DO 260 J = L, M 916 S = 0.0d0 917C 918 DO 250 K = L, N 919 250 S = S + U(J,K) * U(I,K) 920C 921 DO 260 K = L, N 922 U(J,K) = U(J,K) + S * RV1(K) 923 260 CONTINUE 924C 925 270 DO 280 K = L, N 926 280 U(I,K) = SCALE * U(I,K) 927C 928 290 S1 = MAX(S1,ABS(W(I))+ABS(RV1(I))) 929 300 CONTINUE 930C .......... ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS .......... 931 IF (.NOT. MATV) GO TO 410 932C .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... 933 DO 400 II = 1, N 934 I = N + 1 - II 935 IF (I .EQ. N) GO TO 390 936 IF (G .EQ. 0.0d0) GO TO 360 937C 938 DO 320 J = L, N 939C .......... DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW .......... 940 320 V(J,I) = (U(I,J) / U(I,L)) / G 941C 942 DO 350 J = L, N 943 S = 0.0d0 944C 945 DO 340 K = L, N 946 340 S = S + U(I,K) * V(K,J) 947C 948 DO 350 K = L, N 949 V(K,J) = V(K,J) + S * V(K,I) 950 350 CONTINUE 951C 952 360 DO 380 J = L, N 953 V(I,J) = 0.0d0 954 V(J,I) = 0.0d0 955 380 CONTINUE 956C 957 390 V(I,I) = 1.0d0 958 G = RV1(I) 959 L = I 960 400 CONTINUE 961C .......... ACCUMULATION OF LEFT-HAND TRANSFORMATIONS .......... 962 410 IF (.NOT. MATU) GO TO 510 963C ..........FOR I=MIN(M,N) STEP -1 UNTIL 1 DO -- .......... 964 MN = N 965 IF (M .LT. N) MN = M 966C 967 DO 500 II = 1, MN 968 I = MN + 1 - II 969 L = I + 1 970 G = W(I) 971 IF (I .EQ. N) GO TO 430 972C 973 DO 420 J = L, N 974 420 U(I,J) = 0.0d0 975C 976 430 IF (G .EQ. 0.0d0) GO TO 475 977 IF (I .EQ. MN) GO TO 460 978C 979 DO 450 J = L, N 980 S = 0.0d0 981C 982 DO 440 K = L, M 983 440 S = S + U(K,I) * U(K,J) 984C .......... DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW .......... 985 F = (S / U(I,I)) / G 986C 987 DO 450 K = I, M 988 U(K,J) = U(K,J) + F * U(K,I) 989 450 CONTINUE 990C 991 460 DO 470 J = I, M 992 470 U(J,I) = U(J,I) / G 993C 994 GO TO 490 995C 996 475 DO 480 J = I, M 997 480 U(J,I) = 0.0d0 998C 999 490 U(I,I) = U(I,I) + 1.0d0 1000 500 CONTINUE 1001C .......... DIAGONALIZATION OF THE BIDIAGONAL FORM .......... 1002 510 CONTINUE 1003C .......... FOR K=N STEP -1 UNTIL 1 DO -- .......... 1004 DO 700 KK = 1, N 1005 K1 = N - KK 1006 K = K1 + 1 1007 ITS = 0 1008C .......... TEST FOR SPLITTING. 1009C FOR L=K STEP -1 UNTIL 1 DO -- .......... 1010 520 DO 530 LL = 1, K 1011 L1 = K - LL 1012 L = L1 + 1 1013 IF (S1 + ABS(RV1(L)) .EQ. S1) GO TO 565 1014C .......... RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT 1015C THROUGH THE BOTTOM OF THE LOOP .......... 1016 IF (S1 + ABS(W(L1)) .EQ. S1) GO TO 540 1017 530 CONTINUE 1018C .......... CANCELLATION OF RV1(L) IF L GREATER THAN 1 .......... 1019 540 C = 0.0d0 1020 S = 1.0d0 1021C 1022 DO 560 I = L, K 1023 F = S * RV1(I) 1024 RV1(I) = C * RV1(I) 1025 IF (S1 + ABS(F) .EQ. S1) GO TO 565 1026 G = W(I) 1027 H = dsqrt(f*f+g*g) 1028 W(I) = H 1029 C = G / H 1030 S = -F / H 1031 IF (.NOT. MATU) GO TO 560 1032C 1033 DO 550 J = 1, M 1034 Y = U(J,L1) 1035 Z = U(J,I) 1036 U(J,L1) = Y * C + Z * S 1037 U(J,I) = -Y * S + Z * C 1038 550 CONTINUE 1039C 1040 560 CONTINUE 1041C .......... TEST FOR CONVERGENCE .......... 1042 565 Z = W(K) 1043 IF (L .EQ. K) GO TO 650 1044C .......... SHIFT FROM BOTTOM 2 BY 2 MINOR .......... 1045 IF (ITS .EQ. 30) GO TO 1000 1046 ITS = ITS + 1 1047 X = W(L) 1048 Y = W(K1) 1049 G = RV1(K1) 1050 H = RV1(K) 1051 F = 0.5d0 * (((G + Z) / H) * ((G - Z) / Y) + Y / H - H / Y) 1052 G = dsqrt(F*F+1.0d0) 1053 F = X - (Z / X) * Z + (H / X) * (Y / (F + SIGN(G,F)) - H) 1054C .......... NEXT QR TRANSFORMATION .......... 1055 C = 1.0d0 1056 S = 1.0d0 1057C 1058 DO 600 I1 = L, K1 1059 I = I1 + 1 1060 G = RV1(I) 1061 Y = W(I) 1062 H = S * G 1063 G = C * G 1064 Z = dsqrt(F*f+H*h) 1065 RV1(I1) = Z 1066 C = F / Z 1067 S = H / Z 1068 F = X * C + G * S 1069 G = -X * S + G * C 1070 H = Y * S 1071 Y = Y * C 1072 IF (.NOT. MATV) GO TO 575 1073C 1074 DO 570 J = 1, N 1075 X = V(J,I1) 1076 Z = V(J,I) 1077 V(J,I1) = X * C + Z * S 1078 V(J,I) = -X * S + Z * C 1079 570 CONTINUE 1080C 1081 575 Z = dsqrt(F*f+H*h) 1082 W(I1) = Z 1083C .......... ROTATION CAN BE ARBITRARY IF Z IS ZERO .......... 1084 IF (Z .EQ. 0.0d0) GO TO 580 1085 C = F / Z 1086 S = H / Z 1087 580 F = C * G + S * Y 1088 X = -S * G + C * Y 1089 IF (.NOT. MATU) GO TO 600 1090C 1091 DO 590 J = 1, M 1092 Y = U(J,I1) 1093 Z = U(J,I) 1094 U(J,I1) = Y * C + Z * S 1095 U(J,I) = -Y * S + Z * C 1096 590 CONTINUE 1097C 1098 600 CONTINUE 1099C 1100 RV1(L) = 0.0d0 1101 RV1(K) = F 1102 W(K) = X 1103 GO TO 520 1104C .......... CONVERGENCE .......... 1105 650 IF (Z .GE. 0.0d0) GO TO 700 1106C .......... W(K) IS MADE NON-NEGATIVE .......... 1107 W(K) = -Z 1108 IF (.NOT. MATV) GO TO 700 1109C 1110 DO 690 J = 1, N 1111 690 V(J,K) = -V(J,K) 1112C 1113 700 CONTINUE 1114C 1115 GO TO 1001 1116C .......... SET ERROR -- NO CONVERGENCE TO A 1117C SINGULAR VALUE AFTER 30 ITERATIONS .......... 1118 1000 IERR = K 1119 1001 RETURN 1120 END 1121