Millepede-II  V04-00-00_preview
 All Classes Files Functions Variables Enumerator Pages
mpnum.f90
Go to the documentation of this file.
1 
2 ! Code converted using TO_F90 by Alan Miller
3 ! Date: 2012-03-16 Time: 11:07:48
4 
51 
52 !----------------------------------------------------------------------
71 
72 SUBROUTINE sqminv(v,b,n,nrank,diag,next) ! matrix inversion
73  IMPLICIT NONE
74  INTEGER :: i
75  INTEGER :: ij
76  INTEGER :: j
77  INTEGER :: jj
78  INTEGER :: jk
79  INTEGER :: jl
80  INTEGER :: k
81  INTEGER :: kk
82  INTEGER :: l
83  INTEGER :: last
84  INTEGER :: lk
85  INTEGER :: next0
86 
87  DOUBLE PRECISION, INTENT(IN OUT) :: v(*)
88  DOUBLE PRECISION, INTENT(OUT) :: b(n)
89  INTEGER, INTENT(IN) :: n
90  INTEGER, INTENT(OUT) :: nrank
91  DOUBLE PRECISION, INTENT(OUT) :: diag(n)
92  INTEGER, INTENT(OUT) :: next(n)
93  DOUBLE PRECISION :: vkk
94  DOUBLE PRECISION :: vjk
95 
96  DOUBLE PRECISION, PARAMETER :: eps=1.0d-10
97  ! ...
98  next0=1
99  l=1
100  DO i=1,n
101  next(i)=i+1 ! set "next" pointer
102  diag(i)=abs(v((i*i+i)/2)) ! save abs of diagonal elements
103  END DO
104  next(n)=-1 ! end flag
105 
106  nrank=0
107  DO i=1,n ! start of loop
108  k =0
109  vkk=0.0d0
110 
111  j=next0
112  last=0
113 05 IF(j > 0) THEN
114  jj=(j*j+j)/2
115  IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
116  vkk=v(jj)
117  k=j
118  l=last
119  END IF
120  last=j
121  j=next(last)
122  go to 05
123  END IF
124 
125  IF(k /= 0) THEN ! pivot found
126  kk=(k*k+k)/2
127  IF(l == 0) THEN
128  next0=next(k)
129  ELSE
130  next(l)=next(k)
131  END IF
132  next(k)=0 ! index is used, reset
133  nrank=nrank+1 ! increase rank and ...
134  vkk =1.0/vkk
135  v(kk) =-vkk
136  b(k) =b(k)*vkk
137  jk =kk-k
138  jl =0
139  DO j=1,n ! elimination
140  IF(j == k) THEN
141  jk=kk
142  jl=jl+j
143  ELSE
144  IF(j < k) THEN
145  jk=jk+1
146  ELSE
147  jk=jk+j-1
148  END IF
149  vjk =v(jk)
150  v(jk)=vkk*vjk
151  b(j) =b(j)-b(k)*vjk
152  lk =kk-k
153  DO l=1,j
154  jl=jl+1
155  IF(l == k) THEN
156  lk=kk
157  ELSE
158  IF(l < k) THEN
159  lk=lk+1
160  ELSE
161  lk=lk+l-1
162  END IF
163  v(jl)=v(jl)-v(lk)*vjk
164  END IF
165  END DO
166  END IF
167  END DO
168  ELSE
169  DO k=1,n
170  IF(next(k) /= 0) THEN
171  b(k)=0.0d0 ! clear vector element
172  DO j=1,k
173  IF(next(j) /= 0) v((k*k-k)/2+j)=0.0d0 ! clear matrix row/col
174  END DO
175  END IF
176  END DO
177  go to 10
178  END IF
179  END DO ! end of loop
180  10 DO ij=1,(n*n+n)/2
181  v(ij)=-v(ij) ! finally reverse sign of all matrix elements
182  END DO
183 END SUBROUTINE sqminv
184 
197 
198 SUBROUTINE sqminl(v,b,n,nrank,diag,next) !
199  USE mpdef
200 
201  IMPLICIT NONE
202  INTEGER :: i
203  INTEGER :: j
204  INTEGER :: k
205  INTEGER :: l
206  INTEGER :: last
207  INTEGER :: next0
208 
209  DOUBLE PRECISION, INTENT(IN OUT) :: v(*)
210  DOUBLE PRECISION, INTENT(OUT) :: b(n)
211  INTEGER, INTENT(IN) :: n
212  INTEGER, INTENT(OUT) :: nrank
213  DOUBLE PRECISION, INTENT(OUT) :: diag(n)
214  INTEGER, INTENT(OUT) :: next(n)
215  INTEGER*8 :: i8
216  INTEGER*8 :: j8
217  INTEGER*8 :: jj
218  INTEGER*8 :: k8
219  INTEGER*8 :: kk
220  INTEGER*8 :: kkmk
221  INTEGER*8 :: jk
222  INTEGER*8 :: jl
223  INTEGER*8 :: llk
224  INTEGER*8 :: ljl
225 
226  DOUBLE PRECISION :: vkk
227  DOUBLE PRECISION :: vjk
228 
229  DOUBLE PRECISION, PARAMETER :: eps=1.0d-10
230  ! ...
231  next0=1
232  l=1
233  DO i=1,n
234  i8=int8(i)
235  next(i)=i+1 ! set "next" pointer
236  diag(i)=abs(v((i8*i8+i8)/2)) ! save abs of diagonal elements
237  END DO
238  next(n)=-1 ! end flag
239 
240  nrank=0
241  DO i=1,n ! start of loop
242  k =0
243  vkk=0.0d0
244  j=next0
245  last=0
246 05 IF(j > 0) THEN
247  j8=int8(j)
248  jj=(j8*j8+j8)/2
249  IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
250  vkk=v(jj)
251  k=j
252  l=last
253  END IF
254  last=j
255  j=next(last)
256  go to 05
257  END IF
258 
259  IF(k /= 0) THEN ! pivot found
260  k8=int8(k)
261  kk=(k8*k8+k8)/2
262  kkmk=kk-k8
263  IF(l == 0) THEN
264  next0=next(k)
265  ELSE
266  next(l)=next(k)
267  END IF
268  next(k)=0 ! index is used, reset
269  nrank=nrank+1 ! increase rank and ...
270  vkk =1.0/vkk
271  v(kk) =-vkk
272  b(k) =b(k)*vkk
273  ! elimination
274  jk =kkmk
275  DO j=1,n
276  IF(j == k) THEN
277  jk=kk
278  ELSE
279  IF(j < k) THEN
280  jk=jk+1
281  ELSE
282  jk=jk+int8(j)-1
283  END IF
284  v(jk)=v(jk)*vkk
285  END IF
286  END DO
287  ! parallelize row loop
288  ! slot of 128 'J' for next idle thread
289  !$OMP PARALLEL DO &
290  !$OMP PRIVATE(JL,JK,L,LJL,LLK,VJK,J8) &
291  !$OMP SCHEDULE(DYNAMIC,128)
292  DO j=n,1,-1
293  j8=int8(j)
294  jl=j8*(j8-1)/2
295  IF(j /= k) THEN
296  IF(j < k) THEN
297  jk=kkmk+j8
298  ELSE
299  jk=k8+jl
300  END IF
301  vjk =v(jk)/vkk
302  b(j) =b(j)-b(k)*vjk
303  ljl=jl
304  llk=kkmk
305  DO l=1,min(j,k-1)
306  ljl=ljl+1
307  llk=llk+1
308  v(ljl)=v(ljl)-v(llk)*vjk
309  END DO
310  ljl=ljl+1
311  llk=kk
312  DO l=k+1,j
313  ljl=ljl+1
314  llk=llk+l-1
315  v(ljl)=v(ljl)-v(llk)*vjk
316  END DO
317  END IF
318  END DO
319  !$OMP END PARALLEL DO
320  ELSE
321  DO k=1,n
322  k8=int8(k)
323  kk=(k8*k8-k8)/2
324  IF(next(k) /= 0) THEN
325  b(k)=0.0d0 ! clear vector element
326  DO j=1,k
327  IF(next(j) /= 0) v(kk+int8(j))=0.0d0 ! clear matrix row/col
328  END DO
329  END IF
330  END DO
331  go to 10
332  END IF
333  END DO ! end of loop
334  10 DO jj=1,(int8(n)*int8(n)+int8(n))/2
335  v(jj)=-v(jj) ! finally reverse sign of all matrix elements
336  END DO
337 END SUBROUTINE sqminl
338 
350 
351 SUBROUTINE devrot(n,diag,u,v,work,iwork) ! diagonalization
352  IMPLICIT NONE
353 
354  INTEGER, INTENT(IN) :: n
355  DOUBLE PRECISION, INTENT(OUT) :: diag(n)
356  DOUBLE PRECISION, INTENT(OUT) :: u(n,n)
357  DOUBLE PRECISION, INTENT(IN) :: v(*)
358  DOUBLE PRECISION, INTENT(OUT) :: work(n)
359  INTEGER, INTENT(OUT) :: iwork(n)
360 
361 
362  INTEGER, PARAMETER :: itmax=30
363  DOUBLE PRECISION, PARAMETER :: tol=1.0d-16
364  DOUBLE PRECISION, PARAMETER :: eps=1.0d-16
365 
366  DOUBLE PRECISION :: f
367  DOUBLE PRECISION :: g
368  DOUBLE PRECISION :: h
369  DOUBLE PRECISION :: sh
370  DOUBLE PRECISION :: hh
371  DOUBLE PRECISION :: b
372  DOUBLE PRECISION :: p
373  DOUBLE PRECISION :: r
374  DOUBLE PRECISION :: s
375  DOUBLE PRECISION :: c
376  DOUBLE PRECISION :: workd
377 
378  INTEGER :: ij
379  INTEGER :: i
380  INTEGER :: j
381  INTEGER :: k
382  INTEGER :: l
383  INTEGER :: m
384  INTEGER :: ll
385  ! ...
386  ! 1. part: symmetric matrix V reduced to tridiagonal from
387  ij=0
388  DO i=1,n
389  DO j=1,i
390  ij=ij+1
391  u(i,j)=v(ij) ! copy half of symmetric matirx
392  END DO
393  END DO
394 
395  DO i=n,2,-1
396  l=i-2
397  f=u(i,i-1)
398  g=0.0d0
399  IF(l /= 0) THEN
400  DO k=1,l
401  IF(abs(u(i,k)) > tol) g=g+u(i,k)*u(i,k)
402  END DO
403  h=g+f*f
404  END IF
405  IF(g < tol) THEN ! G too small
406  work(i)=f ! skip transformation
407  h =0.0d0
408  ELSE
409  l=l+1
410  sh=sqrt(h)
411  IF(f >= 0.0d0) sh=-sh
412  g=sh
413  work(i)=sh
414  h=h-f*g
415  u(i,i-1)=f-g
416  f=0.0d0
417  DO j=1,l
418  u(j,i)=u(i,j)/h
419  g=0.0d0
420  ! form element of a u
421  DO k=1,j
422  IF(abs(u(j,k)) > tol.AND.abs(u(i,k)) > tol) THEN
423  g=g+u(j,k)*u(i,k)
424  END IF
425  END DO
426  DO k=j+1,l
427  IF(abs(u(k,j)) > tol.AND.abs(u(i,k)) > tol) THEN
428  g=g+u(k,j)*u(i,k)
429  END IF
430  END DO
431  work(j)=g/h
432  f=f+g*u(j,i)
433  END DO
434  ! form k
435  hh=f/(h+h)
436  ! form reduced a
437  DO j=1,l
438  f=u(i,j)
439  work(j)=work(j)-hh*f
440  g=work(j)
441  DO k=1,j
442  u(j,k)=u(j,k)-f*work(k)-g*u(i,k)
443  END DO
444  END DO
445  END IF
446  diag(i)=h
447  END DO
448 
449  diag(1)=0.0d0
450  work(1)=0.0d0
451 
452  ! accumulation of transformation matrices
453  DO i=1,n
454  IF(diag(i) /= 0.0) THEN
455  DO j=1,i-1
456  g=0.0d0
457  DO k=1,i-1
458  g=g+u(i,k)*u(k,j)
459  END DO
460  DO k=1,i-1
461  u(k,j)=u(k,j)-g*u(k,i)
462  END DO
463  END DO
464  END IF
465  diag(i)=u(i,i)
466  u(i,i)=1.0d0
467  DO j=1,i-1
468  u(i,j)=0.0d0
469  u(j,i)=0.0d0
470  END DO
471  END DO
472 
473  ! 2. part: diagonalization of tridiagonal matrix
474  DO i=2,n
475  work(i-1)=work(i)
476  END DO
477  work(n)=0.0d0
478  b=0.0d0
479  f=0.0d0
480 
481  DO l=1,n
482  j=0
483  h=eps*(abs(diag(l))+abs(work(l)))
484  IF(b < h) b=h
485  DO m=l,n
486  IF(abs(work(m)) <= b) go to 10 ! look for small sub-diagonal element
487  END DO
488  m=l
489 10 IF(m == l) go to 30
490  ! next iteration
491 20 IF(j == itmax) THEN
492  WRITE(*,*) 'DEVROT: Iteration limit reached'
493  stop
494  END IF
495  j=j+1
496  g=diag(l)
497  p=(diag(l+1)-g)/(2.0d0*work(l))
498  r=sqrt(1.0d0+p*p)
499  diag(l)=work(l)
500  IF(p < 0.0d0) diag(l)=diag(l)/(p-r)
501  IF(p >= 0.0d0) diag(l)=diag(l)/(p+r)
502  h=g-diag(l)
503  DO i=l+1,n
504  diag(i)=diag(i)-h
505  END DO
506  f=f+h
507  ! QL transformation
508  p=diag(m)
509  c=1.0d0
510  s=0.0d0
511  DO i=m-1,l,-1 ! reverse loop
512  g=c*work(i)
513  h=c*p
514  IF(abs(p) >= abs(work(i))) THEN
515  c=work(i)/p
516  r=sqrt(1.0d0+c*c)
517  work(i+1)=s*p*r
518  s=c/r
519  c=1.0d0/r
520  ELSE
521  c=p/work(i)
522  r=sqrt(1.0d0+c*c)
523  work(i+1)=s*work(i)*r
524  s=1.0d0/r
525  c=c/r
526  END IF
527  p=c*diag(i)-s*g
528  diag(i+1)=h+s*(c*g+s*diag(i))
529  ! form vector
530  DO k=1,n
531  h=u(k,i+1)
532  u(k,i+1)=s*u(k,i)+c*h
533  u(k,i)=c*u(k,i)-s*h
534  END DO
535  END DO
536  work(l)=s*p
537  diag(l)=c*p
538  IF(abs(work(l)) > b) go to 20 ! next iteration
539 30 diag(l)=diag(l)+f
540  END DO
541  DO i=1,n
542  iwork(i)=i
543  END DO
544 
545  m=1
546 40 m=1+3*m ! determine initial increment
547  IF(m <= n) go to 40
548 50 m=m/3
549  DO j=1,n-m ! sort with increment M
550  l=j
551 60 IF(diag(iwork(l+m)) > diag(iwork(l))) THEN ! compare
552  ll=iwork(l+m) ! exchange the two index values
553  iwork(l+m)=iwork(l)
554  iwork(l)=ll
555  l=l-m
556  IF(l > 0) go to 60
557  END IF
558  END DO
559  IF(m > 1) go to 50
560 
561  DO i=1,n
562  IF(iwork(i) /= i) THEN
563  ! move vector from position I to the work area
564  workd=diag(i)
565  DO l=1,n
566  work(l)=u(l,i)
567  END DO
568  k=i
569 70 j=k
570  k=iwork(j)
571  iwork(j)=j
572  IF(k /= i) THEN
573  ! move vector from position K to the (free) position J
574  diag(j)=diag(k)
575  DO l=1,n
576  u(l,j)=u(l,k)
577  END DO
578  go to 70
579  END IF
580  ! move vector from the work area to position J
581  diag(j)=workd
582  DO l=1,n
583  u(l,j)=work(l)
584  END DO
585  END IF
586  END DO
587 END SUBROUTINE devrot
588 
590 SUBROUTINE devsig(n,diag,u,b,coef)
591  IMPLICIT NONE
592 
593  INTEGER, INTENT(IN) :: n
594  DOUBLE PRECISION, INTENT(IN) :: diag(n)
595  DOUBLE PRECISION, INTENT(IN) :: u(n,n)
596  DOUBLE PRECISION, INTENT(IN) :: b(n)
597  DOUBLE PRECISION, INTENT(OUT) :: coef(n)
598  INTEGER :: i
599  INTEGER :: j
600  DOUBLE PRECISION :: dsum
601  ! ...
602  DO i=1,n
603  coef(i)=0.0d0
604  IF(diag(i) > 0.0d0) THEN
605  dsum=0.0d0
606  DO j=1,n
607  dsum=dsum+u(j,i)*b(j)
608  END DO
609  coef(i)=abs(dsum)/sqrt(diag(i))
610  END IF
611  END DO
612 END SUBROUTINE devsig
613 
614 
625 
626 SUBROUTINE devsol(n,diag,u,b,x,work)
627  IMPLICIT NONE
628 
629  INTEGER, INTENT(IN) :: n
630  DOUBLE PRECISION, INTENT(IN) :: diag(n)
631  DOUBLE PRECISION, INTENT(IN) :: u(n,n)
632  DOUBLE PRECISION, INTENT(IN) :: b(n)
633  DOUBLE PRECISION, INTENT(OUT) :: x(n)
634  DOUBLE PRECISION, INTENT(OUT) :: work(n)
635  INTEGER :: i
636  INTEGER :: j
637  INTEGER :: jj
638  DOUBLE PRECISION :: s
639  ! ...
640  DO j=1,n
641  s=0.0d0
642  work(j)=0.0d0
643  IF(diag(j) /= 0.0d0) THEN
644  DO i=1,n
645  ! j-th eigenvector is U(.,J)
646  s=s+u(i,j)*b(i)
647  END DO
648  work(j)=s/diag(j)
649  END IF
650  END DO
651 
652  DO j=1,n
653  s=0.0d0
654  DO jj=1,n
655  s=s+u(j,jj)*work(jj)
656  END DO
657  x(j)=s
658  END DO
659 ! WRITE(*,*) 'DEVSOL'
660 ! WRITE(*,*) 'X ',X
661 END SUBROUTINE devsol
662 
670 
671 SUBROUTINE devinv(n,diag,u,v)
672 
673  IMPLICIT NONE
674  INTEGER :: i
675  INTEGER :: ij
676  INTEGER :: j
677  INTEGER :: k
678 
679  INTEGER, INTENT(IN) :: n
680  DOUBLE PRECISION, INTENT(IN) :: diag(n)
681  DOUBLE PRECISION, INTENT(IN) :: u(n,n)
682  DOUBLE PRECISION, INTENT(OUT) :: v(*)
683  DOUBLE PRECISION:: dsum
684  ! ...
685  ij=0
686  DO i=1,n
687  DO j=1,i
688  ij=ij+1
689  dsum=0.0d0
690  DO k=1,n
691  IF(diag(k) /= 0.0d0) THEN
692  dsum=dsum+u(i,k)*u(j,k)/diag(k)
693  END IF
694  END DO
695  v(ij)=dsum
696  END DO
697  END DO
698 END SUBROUTINE devinv
699 
700 
719 SUBROUTINE choldc(g,n)
720  IMPLICIT NONE
721  INTEGER :: i
722  INTEGER :: ii
723  INTEGER :: j
724  INTEGER :: jj
725  INTEGER :: k
726  INTEGER :: kk
727 
728  DOUBLE PRECISION, INTENT(IN OUT) :: g(*)
729  INTEGER, INTENT(IN) :: n
730 
731  DOUBLE PRECISION :: ratio
732  ! ...
733  ii=0
734  DO i=1,n
735  ii=ii+i
736  IF(g(ii) /= 0.0) g(ii)=1.0/g(ii) ! (I,I) div !
737  jj=ii
738  DO j=i+1,n
739  ratio=g(i+jj)*g(ii) ! (I,J) (I,I)
740  kk=jj
741  DO k=j,n
742  g(kk+j)=g(kk+j)-g(kk+i)*ratio ! (K,J) (K,I)
743  kk=kk+k
744  END DO ! K
745  g(i+jj)=ratio ! (I,J)
746  jj=jj+j
747  END DO ! J
748  END DO ! I
749  return
750 END SUBROUTINE choldc
751 
763 SUBROUTINE cholsl(g,x,n)
764  IMPLICIT NONE
765  DOUBLE PRECISION :: dsum
766  INTEGER :: i
767  INTEGER :: ii
768  INTEGER :: k
769  INTEGER :: kk
770 
771  DOUBLE PRECISION, INTENT(IN) :: g(*)
772  DOUBLE PRECISION, INTENT(IN OUT) :: x(n)
773  INTEGER, INTENT(IN) :: n
774 
775  ii=0
776  DO i=1,n
777  dsum=x(i)
778  DO k=1,i-1
779  dsum=dsum-g(k+ii)*x(k) ! (K,I)
780  END DO
781  x(i)=dsum
782  ii=ii+i
783  END DO
784  DO i=n,1,-1
785  dsum=x(i)*g(ii) ! (I,I)
786  kk=ii
787  DO k=i+1,n
788  dsum=dsum-g(kk+i)*x(k) ! (K,I)
789  kk=kk+k
790  END DO
791  x(i)=dsum
792  ii=ii-i
793  END DO
794  return
795 END SUBROUTINE cholsl
796 
807 SUBROUTINE cholin(g,v,n)
808  IMPLICIT NONE
809  DOUBLE PRECISION :: dsum
810  INTEGER :: i
811  INTEGER :: ii
812  INTEGER :: j
813  INTEGER :: k
814  INTEGER :: l
815  INTEGER :: m
816 
817  DOUBLE PRECISION, INTENT(IN) :: g(*)
818  DOUBLE PRECISION, INTENT( OUT) :: v(*)
819  INTEGER, INTENT(IN) :: n
820 
821  ii=(n*n-n)/2
822  DO i=n,1,-1
823  dsum=g(ii+i) ! (I,I)
824  DO j=i,1,-1
825  DO k=j+1,n
826  l=min(i,k)
827  m=max(i,k)
828  dsum=dsum-g(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K)
829  END DO
830  v(ii+j)=dsum ! (I,J)
831  dsum=0.0d0
832  END DO
833  ii=ii-i+1
834  END DO
835 END SUBROUTINE cholin
836 
837 ! variable band matrix operations ----------------------------------
838 
859 
860 SUBROUTINE vabdec(n,val,ilptr)
861  IMPLICIT NONE
862 
863  INTEGER :: i
864  INTEGER :: in
865  INTEGER :: j
866  INTEGER :: k
867  INTEGER :: kj
868  INTEGER :: mj
869  INTEGER :: mk
870  DOUBLE PRECISION :: sn
871  DOUBLE PRECISION :: beta
872  DOUBLE PRECISION :: delta
873  DOUBLE PRECISION :: theta
874 
875  INTEGER, INTENT(IN) :: n
876  DOUBLE PRECISION, INTENT(IN OUT) :: val(*)
877  INTEGER, INTENT(IN) :: ilptr(n)
878 
879  DOUBLE PRECISION :: dgamma
880  DOUBLE PRECISION :: xi
881  DOUBLE PRECISION :: eps
882  DOUBLE PRECISION :: prd
883  DOUBLE PRECISION :: valkj
884 
885  DOUBLE PRECISION, PARAMETER :: one=1.0d0
886  DOUBLE PRECISION, PARAMETER :: two=2.0d0
887  ! DATA EPS/0.0D0/
888  DATA eps/2.22044605e-16/
889  SAVE eps
890  ! ...
891  IF(eps == 0.0d0) THEN
892  eps = two**(-12)
893 10 eps = eps/two
894  prd=one
895  prd=prd+one*eps ! CALL dbaxpy(1,one,eps,prd)
896  IF(prd > one) go to 10
897  eps=eps*two ! EPS is machine presision
898  WRITE(*,*) 'Machine precision is ',eps
899  END IF
900 
901  WRITE(*,*) 'Variable band matrix Cholesky decomposition'
902 
903  dgamma=0.0d0
904  i=1
905  DO j=1,ilptr(n) ! loop thrugh all matrix elements
906  IF(ilptr(i) == j) THEN ! diagonal element
907  IF(val(j) <= 0.0d0) go to 01 ! exit loop for negative diag
908  dgamma=max(dgamma,abs(val(j))) ! max diagonal element
909  i=i+1
910  END IF
911  END DO
912  i=n+1
913 01 in=i-1 ! IN positive diagonal elements
914  WRITE(*,*) ' ',in,' positive diagonal elements'
915  xi=0.0d0
916  i=1
917  DO j=1,ilptr(in) ! loop for positive diagonal elements
918  ! through all matrix elements
919  IF(ilptr(i) == j) THEN ! diagonal element
920  i=i+1
921  ELSE
922  xi=max(xi,abs(val(j))) ! Xi = abs(max) off-diagonal element
923  END IF
924  END DO
925 
926  delta=eps*max(1.0d0,dgamma+xi)
927  sn=1.0d0
928  IF(n > 1) sn=1.0d0/sqrt(dfloat(n*n-1))
929  beta=sqrt(max(eps,dgamma,xi*sn)) ! beta
930  WRITE(*,*) ' DELTA and BETA ',delta,beta
931 
932  DO k=2,n
933  mk=k-ilptr(k)+ilptr(k-1)+1
934 
935  theta=0.0d0
936 
937  DO j=mk,k
938  mj=j-ilptr(j)+ilptr(j-1)+1
939  kj=ilptr(k)-k+j ! index kj
940 
941  DO i=max(mj,mk),j-1
942  val(kj)=val(kj) & ! L_kj := L_kj - L_ki D_ii L_ji
943  -val(ilptr(k)-k+i)*val(ilptr(i))*val(ilptr(j)-j+i)
944 
945  END DO !
946 
947  theta=max(theta,abs(val(kj))) ! maximum value of row
948 
949  IF(j /= k) THEN
950  IF(val(ilptr(j)) /= 0.0d0) THEN
951  val(kj)=val(kj)/val(ilptr(j))
952  ELSE
953  val(kj)=0.0d0
954  END IF
955  END IF ! L_kj := L_kj/D_jj ! D_kk
956 
957  IF(j == k) THEN
958  valkj=val(kj)
959  IF(k <= in) THEN
960  val(kj)=max(abs(val(kj)),(theta/beta)**2,delta)
961  IF(valkj /= val(kj)) THEN
962  WRITE(*,*) ' Index K=',k
963  WRITE(*,*) ' ',valkj,val(kj), (theta/beta)**2,delta,theta
964  END IF
965  END IF
966  END IF
967  END DO ! J
968 
969  END DO ! K
970 
971  DO k=1,n
972  IF(val(ilptr(k)) /= 0.0d0) val(ilptr(k))=1.0d0/val(ilptr(k))
973  END DO
974 
975  return
976 END SUBROUTINE vabdec
977 
983 
984 SUBROUTINE vabmmm(n,val,ilptr)
985  IMPLICIT NONE
986  INTEGER :: k
987  INTEGER :: kp
988  INTEGER :: kr
989  INTEGER :: ks
990  INTEGER, INTENT(IN) :: n
991  DOUBLE PRECISION, INTENT(IN OUT) :: val(*)
992  INTEGER, INTENT(IN) :: ilptr(n)
993  kr=1
994  ks=1
995  kp=1
996  DO k=1,n
997  IF(val(ilptr(k)) > val(ilptr(ks))) ks=k
998  IF(val(ilptr(k)) < val(ilptr(kr))) kr=k
999  IF(val(ilptr(k)) > 0.0.AND.val(ilptr(k)) < val(ilptr(kp))) kp=k
1000  END DO
1001  WRITE(*,*) ' Index value ',ks,val(ilptr(ks))
1002  WRITE(*,*) ' Index value ',kp,val(ilptr(kp))
1003  WRITE(*,*) ' Index value ',kr,val(ilptr(kr))
1004 
1005  return
1006 END SUBROUTINE vabmmm
1007 
1018 
1019 SUBROUTINE vabslv(n,val,ilptr,x)
1020  IMPLICIT NONE
1021  INTEGER :: j
1022  INTEGER :: k
1023  INTEGER :: mk
1024 
1025  INTEGER, INTENT(IN) :: n
1026  DOUBLE PRECISION, INTENT(IN OUT) :: val(*)
1027  INTEGER, INTENT(IN) :: ilptr(n)
1028  DOUBLE PRECISION, INTENT(IN OUT) :: x(n)
1029  ! ...
1030  DO k=1,n ! forward loop
1031  mk=k-ilptr(k)+ilptr(k-1)+1
1032  DO j=mk,k-1
1033  x(k)=x(k)-val(ilptr(k)-k+j)*x(j) ! X_k := X_k - L_kj B_j
1034  END DO
1035  END DO ! K
1036 
1037  DO k=1,n ! divide by diagonal elements
1038  x(k)=x(k)*val(ilptr(k)) ! X_k := X_k*D_kk
1039  END DO
1040 
1041  DO k=n,1,-1 ! backward loop
1042  mk=k-ilptr(k)+ilptr(k-1)+1
1043  DO j=mk,k-1
1044  x(j)=x(j)-val(ilptr(k)-k+j)*x(k) ! X_j := X_j - L_kj X_k
1045  END DO
1046  END DO ! K
1047 END SUBROUTINE vabslv
1048 
1049 ! matrix/vector products -------------------------------------------
1050 
1057 
1058 DOUBLE PRECISION FUNCTION dbdot(n,dx,dy)
1059  IMPLICIT NONE
1060  INTEGER:: i
1061  DOUBLE PRECISION :: dtemp
1062 
1063  INTEGER, INTENT(IN) :: n
1064  DOUBLE PRECISION, INTENT(IN) :: dx(*)
1065  DOUBLE PRECISION, INTENT(IN) :: dy(*)
1066  ! ...
1067  dtemp=0.0d0
1068  DO i = 1,mod(n,5)
1069  dtemp=dtemp+dx(i)*dy(i)
1070  END DO
1071  DO i =mod(n,5)+1,n,5
1072  dtemp=dtemp+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2) &
1073  +dx(i+3)*dy(i+3)+dx(i+4)*dy(i+4)
1074  END DO
1075  dbdot=dtemp
1076 END FUNCTION dbdot
1077 
1086 
1087 SUBROUTINE dbaxpy(n,da,dx,dy)
1088  IMPLICIT NONE
1089  INTEGER:: i
1090 
1091  INTEGER, INTENT(IN) :: n
1092  DOUBLE PRECISION, INTENT(IN) :: dx(*)
1093  DOUBLE PRECISION, INTENT(IN OUT) :: dy(*)
1094  DOUBLE PRECISION, INTENT(IN) :: da
1095  ! ...
1096  DO i=1,mod(n,4)
1097  dy(i)=dy(i)+da*dx(i)
1098  END DO
1099  DO i=mod(n,4)+1,n,4
1100  dy(i )=dy(i )+da*dx(i )
1101  dy(i+1)=dy(i+1)+da*dx(i+1)
1102  dy(i+2)=dy(i+2)+da*dx(i+2)
1103  dy(i+3)=dy(i+3)+da*dx(i+3)
1104  END DO
1105 END SUBROUTINE dbaxpy
1106 
1115 
1116 SUBROUTINE dbsvx(v,a,b,n) !
1117  IMPLICIT NONE
1118  INTEGER:: i
1119  INTEGER:: ij
1120  INTEGER:: ijs
1121  INTEGER:: j
1122 
1123  ! B := V * A
1124  ! N N*N N
1125 
1126  INTEGER, INTENT(IN) :: n
1127  DOUBLE PRECISION, INTENT(IN) :: v(*)
1128  DOUBLE PRECISION, INTENT(IN) :: a(*)
1129  DOUBLE PRECISION, INTENT(OUT) :: b(*)
1130 
1131  DOUBLE PRECISION:: dsum
1132  ijs=1
1133  DO i=1,n
1134  dsum=0.0
1135  ij=ijs
1136  DO j=1,n
1137  dsum=dsum+v(ij)*a(j)
1138  IF(j < i) THEN
1139  ij=ij+1
1140  ELSE
1141  ij=ij+j
1142  END IF
1143  END DO
1144  b(i)=dsum
1145  ijs=ijs+i
1146  END DO
1147 END SUBROUTINE dbsvx
1148 
1157 
1158 SUBROUTINE dbsvxl(v,a,b,n) ! LARGE symm. matrix, vector
1159  IMPLICIT NONE
1160  INTEGER:: i
1161  INTEGER:: j
1162 
1163  ! B := V * A
1164  ! N N*N N
1165 
1166  INTEGER, INTENT(IN) :: n
1167  DOUBLE PRECISION, INTENT(IN) :: v(*)
1168  DOUBLE PRECISION, INTENT(IN) :: a(*)
1169  DOUBLE PRECISION, INTENT(OUT) :: b(*)
1170 
1171  DOUBLE PRECISION:: dsum
1172  INTEGER*8 :: ij
1173  INTEGER*8 :: ijs
1174  ijs=1
1175  DO i=1,n
1176  dsum=0.0
1177  ij=ijs
1178  DO j=1,n
1179  dsum=dsum+v(ij)*a(j)
1180  IF(j < i) THEN
1181  ij=ij+1
1182  ELSE
1183  ij=ij+int8(j)
1184  END IF
1185  END DO
1186  b(i)=dsum
1187  ijs=ijs+int8(i)
1188  END DO
1189 END SUBROUTINE dbsvxl
1190 
1198 
1199 SUBROUTINE dbgax(a,x,y,m,n)
1200  IMPLICIT NONE
1201  INTEGER :: i
1202  INTEGER :: ij
1203  INTEGER :: j
1204 
1205  DOUBLE PRECISION, INTENT(IN) :: a(*)
1206  DOUBLE PRECISION, INTENT(IN) :: x(*)
1207  DOUBLE PRECISION, INTENT(OUT) :: y(*)
1208  INTEGER, INTENT(IN) :: m
1209  INTEGER, INTENT(IN) :: n
1210 
1211  ! ...
1212  ij=0
1213  DO i=1,m
1214  y(i)=0.0d0
1215  DO j=1,n
1216  ij=ij+1
1217  y(i)=y(i)+a(ij)*x(j)
1218  END DO
1219  END DO
1220 END SUBROUTINE dbgax
1221 
1234 SUBROUTINE dbavat(v,a,w,n,ms)
1235  IMPLICIT NONE
1236  INTEGER :: i
1237  INTEGER :: ij
1238  INTEGER :: ijs
1239  INTEGER :: il
1240  INTEGER :: j
1241  INTEGER :: jk
1242  INTEGER :: k
1243  INTEGER :: l
1244  INTEGER :: lk
1245  INTEGER :: lkl
1246  INTEGER :: m
1247 
1248  DOUBLE PRECISION, INTENT(IN) :: v(*)
1249  DOUBLE PRECISION, INTENT(IN) :: a(*)
1250  DOUBLE PRECISION, INTENT(INOUT) :: w(*)
1251  INTEGER, INTENT(IN) :: n
1252  INTEGER, INTENT(IN) :: ms
1253 
1254  DOUBLE PRECISION :: cik
1255  ! ...
1256  m=ms
1257  IF (m > 0) THEN
1258  DO i=1,(m*m+m)/2
1259  w(i)=0.0d0 ! reset output matrix
1260  END DO
1261  ELSE
1262  m=-m
1263  END IF
1264 
1265  il=-n
1266  ijs=0
1267  DO i=1,m ! do I
1268  ijs=ijs+i-1 !
1269  il=il+n !
1270  lkl=0 !
1271  DO k=1,n ! do K
1272  cik=0.0d0 !
1273  lkl=lkl+k-1 !
1274  lk=lkl !
1275  DO l=1,k ! do L
1276  lk=lk+1 ! .
1277  cik=cik+a(il+l)*v(lk) ! .
1278  END DO ! end do L
1279  DO l=k+1,n ! do L
1280  lk=lk+l-1 ! .
1281  cik=cik+a(il+l)*v(lk) ! .
1282  END DO ! end do L
1283  jk=k !
1284  ij=ijs !
1285  DO j=1,i ! do J
1286  ij=ij+1 ! .
1287  w(ij)=w(ij)+cik*a(jk) ! .
1288  jk=jk+n ! .
1289  END DO ! end do J
1290  END DO ! end do K
1291  END DO ! end do I
1292 END SUBROUTINE dbavat
1293 
1303 
1304 SUBROUTINE dbmprv(lun,x,v,n)
1305  IMPLICIT NONE
1306  INTEGER :: i
1307  INTEGER :: ii
1308  INTEGER :: ij
1309  INTEGER :: j
1310  INTEGER :: jj
1311  INTEGER :: l
1312  INTEGER :: m
1313  INTEGER :: mc(15)
1314  REAL :: pd
1315  REAL :: rho
1316  REAL :: err
1317 
1318  INTEGER, INTENT(IN) :: lun
1319  DOUBLE PRECISION, INTENT(IN) :: x(*)
1320  DOUBLE PRECISION, INTENT(IN) :: v(*)
1321  INTEGER, INTENT(IN) :: n
1322 
1323  WRITE(lun,103)
1324  WRITE(lun,101)
1325  ii=0
1326  DO i=1,n
1327  ij=ii
1328  ii=ii+i
1329  err=0.0
1330  IF(v(ii) > 0.0) err=sqrt(REAL(v(ii)))
1331  l=0
1332  jj=0
1333  DO j=1,i
1334  jj=jj+j
1335  ij=ij+1
1336  rho=0.0
1337  pd=REAL(v(ii)*v(jj))
1338  IF(pd > 0.0) rho=REAL(v(ij))/sqrt(pd)
1339  l=l+1
1340  mc(l)=int(100.0*abs(rho)+0.5)
1341  IF(rho < 0.0) mc(l)=-mc(l)
1342  IF(j == i.OR.l == 15) THEN
1343  IF(j <= 15) THEN
1344  IF(j == i) THEN
1345  WRITE(lun,102) i,x(i),err,(mc(m),m=1,l-1)
1346  ELSE
1347  WRITE(lun,102) i,x(i),err,(mc(m),m=1,l)
1348  END IF
1349  ELSE
1350  IF(j == i) THEN
1351  WRITE(lun,103) (mc(m),m=1,l-1)
1352  ELSE
1353  WRITE(lun,103) (mc(m),m=1,l)
1354  END IF
1355  l=0
1356  END IF
1357  END IF
1358  END DO
1359  END DO
1360  WRITE(lun,104)
1361  ! 100 RETURN
1362  return
1363 101 format(9x,'Param',7x,'error',7x,'correlation coefficients'/)
1364 102 format(1x,i5,2g12.4,1x,15i5)
1365 103 format(31x,15i5)
1366 104 format(33x,'(correlation coefficients in percent)')
1367 END SUBROUTINE dbmprv
1368 
1376 
1377 SUBROUTINE dbprv(lun,v,n)
1378  IMPLICIT NONE
1379  INTEGER, PARAMETER :: istp=6
1380  INTEGER :: i
1381  INTEGER :: ip
1382  INTEGER :: ipe
1383  INTEGER :: ipn
1384  INTEGER :: ips
1385  INTEGER :: k
1386 
1387  INTEGER, INTENT(IN) :: lun
1388  DOUBLE PRECISION, INTENT(IN) :: v(*)
1389  INTEGER, INTENT(IN) :: n
1390 
1391  WRITE(lun,101)
1392 
1393  DO i=1,n
1394  ips=(i*i-i)/2
1395  ipe=ips+i
1396  ip =ips
1397 100 continue
1398  ipn=ip+istp
1399  WRITE(lun,102), i, ip+1-ips, (v(k),k=ip+1,min(ipn,ipe))
1400  IF (ipn < ipe) THEN
1401  ip=ipn
1402  go to 100
1403  END IF
1404 END DO
1405 return
1406 101 format(1x,'--- DBPRV -----------------------------------')
1407 102 format(1x,2i3,6g12.4)
1408 END SUBROUTINE dbprv
1409 
1410 ! sort -------------------------------------------------------------
1411 
1418 
1419 SUBROUTINE heapf(a,n)
1420  IMPLICIT NONE
1421  !
1422  INTEGER :: i
1423  INTEGER :: j
1424  INTEGER :: l
1425  INTEGER :: r
1426  REAL :: at ! pivot key value
1427 
1428  REAL, INTENT(IN OUT) :: a(*)
1429  INTEGER, INTENT(IN) :: n
1430  ! ...
1431  IF(n <= 1) return
1432  l=n/2+1
1433  r=n
1434 10 IF(l > 1) THEN
1435  l=l-1
1436  at =a(l)
1437  ELSE
1438  at =a(r)
1439  a(r)=a(1)
1440  r=r-1
1441  IF(r == 1) THEN
1442  a(1)=at
1443  return
1444  END IF
1445  END IF
1446  i=l
1447  j=l+l
1448 20 IF(j <= r) THEN
1449  IF(j < r) THEN
1450  IF(a(j) < a(j+1)) j=j+1
1451  END IF
1452  IF(at < a(j)) THEN
1453  a(i)=a(j)
1454  i=j
1455  j=j+j
1456  ELSE
1457  j=r+1
1458  END IF
1459  go to 20
1460  END IF
1461  a(i)=at
1462  go to 10
1463 END SUBROUTINE heapf
1464 
1471 
1472 SUBROUTINE sort1k(a,n)
1473  IMPLICIT NONE
1474  INTEGER :: nlev ! stack size
1475  parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1476  INTEGER :: i
1477  INTEGER :: j
1478  INTEGER :: l
1479  INTEGER :: r
1480  INTEGER :: lev
1481  INTEGER :: lr(nlev)
1482  INTEGER :: lrh
1483  INTEGER :: maxlev
1484  INTEGER :: a1 ! pivot key
1485  INTEGER :: at ! pivot key
1486 
1487  INTEGER, INTENT(IN OUT) :: a(*)
1488  INTEGER, INTENT(IN) :: n
1489  ! ...
1490  IF (n <= 0) return
1491  maxlev=0
1492  lev=0
1493  l=1
1494  r=n
1495 10 IF(r-l == 1) THEN ! sort two elements L and R
1496  IF(a(l) > a(r)) THEN
1497  at=a(l) ! exchange L <-> R
1498  a(l)=a(r)
1499  a(r)=at
1500  END IF
1501  r=l
1502  END IF
1503  IF(r == l) THEN
1504  IF(lev <= 0) THEN
1505  ! WRITE(*,*) 'SORT1K (quicksort): maxlevel used/available =',
1506  ! + MAXLEV,'/64'
1507  return
1508  END IF
1509  lev=lev-2
1510  l=lr(lev+1)
1511  r=lr(lev+2)
1512  ELSE
1513  ! LRH=(L+R)/2
1514  lrh=(l/2)+(r/2) ! avoid bit overflow
1515  IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1516  a1=a(lrh) ! middle
1517  i=l-1 ! find limits [J,I] with [L,R]
1518  j=r+1
1519 20 i=i+1
1520  IF(a(i) < a1) go to 20
1521 30 j=j-1
1522  IF(a(j) > a1) go to 30
1523  IF(i <= j) THEN
1524  at=a(i) ! exchange I <-> J
1525  a(i)=a(j)
1526  a(j)=at
1527  go to 20
1528  END IF
1529  IF(lev+2 > nlev) stop 'SORT1K (quicksort): stack overflow'
1530  IF(r-i < j-l) THEN
1531  lr(lev+1)=l
1532  lr(lev+2)=j
1533  l=i
1534  ELSE
1535  lr(lev+1)=i
1536  lr(lev+2)=r
1537  r=j
1538  END IF
1539  lev=lev+2
1540  maxlev=max(maxlev,lev)
1541  END IF
1542  go to 10
1543 END SUBROUTINE sort1k
1544 
1551 
1552 SUBROUTINE sort2k(a,n)
1553  IMPLICIT NONE
1554  INTEGER:: nlev ! stack size
1555  parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1556  INTEGER:: i
1557  INTEGER::j
1558  INTEGER::l
1559  INTEGER::r
1560  INTEGER::lev
1561  INTEGER::lr(nlev)
1562  INTEGER::lrh
1563  INTEGER::maxlev
1564  INTEGER::a1 ! pivot key
1565  INTEGER::a2 ! pivot key
1566  INTEGER::at ! pivot key
1567 
1568  INTEGER, INTENT(IN OUT) :: a(2,*)
1569  INTEGER, INTENT(IN) :: n
1570  ! ...
1571  maxlev=0
1572  lev=0
1573  l=1
1574  r=n
1575 10 IF(r-l == 1) THEN ! sort two elements L and R
1576  IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
1577  at=a(1,l) ! exchange L <-> R
1578  a(1,l)=a(1,r)
1579  a(1,r)=at
1580  at=a(2,l)
1581  a(2,l)=a(2,r)
1582  a(2,r)=at
1583  END IF
1584  r=l
1585  END IF
1586  IF(r == l) THEN
1587  IF(lev <= 0) THEN
1588  WRITE(*,*) 'SORT2K (quicksort): maxlevel used/available =', maxlev,'/64'
1589  return
1590  END IF
1591  lev=lev-2
1592  l=lr(lev+1)
1593  r=lr(lev+2)
1594  ELSE
1595  ! LRH=(L+R)/2
1596  lrh=(l/2)+(r/2) ! avoid bit overflow
1597  IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1598  a1=a(1,lrh) ! middle
1599  a2=a(2,lrh)
1600  i=l-1 ! find limits [J,I] with [L,R]
1601  j=r+1
1602 20 i=i+1
1603  IF(a(1,i) < a1) go to 20
1604  IF(a(1,i) == a1.AND.a(2,i) < a2) go to 20
1605 30 j=j-1
1606  IF(a(1,j) > a1) go to 30
1607  IF(a(1,j) == a1.AND.a(2,j) > a2) go to 30
1608  IF(i <= j) THEN
1609  at=a(1,i) ! exchange I <-> J
1610  a(1,i)=a(1,j)
1611  a(1,j)=at
1612  at=a(2,i)
1613  a(2,i)=a(2,j)
1614  a(2,j)=at
1615  go to 20
1616  END IF
1617  IF(lev+2 > nlev) stop 'SORT2K (quicksort): stack overflow'
1618  IF(r-i < j-l) THEN
1619  lr(lev+1)=l
1620  lr(lev+2)=j
1621  l=i
1622  ELSE
1623  lr(lev+1)=i
1624  lr(lev+2)=r
1625  r=j
1626  END IF
1627  lev=lev+2
1628  maxlev=max(maxlev,lev)
1629  END IF
1630  go to 10
1631 END SUBROUTINE sort2k
1632 
1640 
1641 REAL FUNCTION chindl(n,nd)
1642  IMPLICIT NONE
1643  INTEGER :: m
1644  INTEGER, INTENT(IN) :: n
1645  INTEGER, INTENT(IN) :: nd
1646  !
1647  REAL:: sn(3)
1648  REAL::table(30,3)
1649  ! REAL PN(3)
1650  ! DATA PN/0.31731,0.0455002785,2.69985E-3/ ! probabilities
1651  DATA sn/0.47523,1.690140,2.782170/
1652  DATA table/ 1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, &
1653  1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, &
1654  1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, &
1655  1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040, &
1656  4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, &
1657  1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, &
1658  1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, &
1659  1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742, &
1660  9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, &
1661  2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, &
1662  2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, &
1663  1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681/
1664  SAVE sn,table
1665  ! ...
1666  IF(nd < 1) THEN
1667  chindl=0.0
1668  ELSE
1669  m=max(1,min(n,3)) ! 1, 2 or 3 sigmas
1670  IF(nd <= 30) THEN
1671  chindl=table(nd,m) ! from table
1672  ELSE ! approximation for ND > 30
1673  chindl=(sn(m)+sqrt(float(nd+nd-1)))**2/float(nd+nd)
1674  END IF
1675  END IF
1676 END FUNCTION chindl
1677 
1711 
1712 SUBROUTINE lltdec(n,c,india,nrkd)
1713  IMPLICIT NONE
1714  INTEGER :: i
1715  INTEGER :: j
1716  INTEGER :: k
1717  INTEGER :: kj
1718  INTEGER :: mj
1719  INTEGER :: mk
1720  DOUBLE PRECISION::diag
1721 
1722  INTEGER, INTENT(IN) :: n
1723  DOUBLE PRECISION, INTENT(IN OUT) :: c(*)
1724  INTEGER, INTENT(IN) :: india(n)
1725  INTEGER, INTENT(OUT) :: nrkd
1726 
1727  ! ...
1728  nrkd=0
1729  diag=0.0d0
1730  IF(c(india(1)) > 0.0) THEN
1731  c(india(1))=1.0d0/sqrt(c(india(1))) ! square root
1732  ELSE
1733  c(india(1))=0.0d0
1734  nrkd=-1
1735  END IF
1736 
1737  DO k=2,n
1738  mk=k-india(k)+india(k-1)+1 ! first index in row K
1739  DO j=mk,k ! loop over row K with index J
1740  mj=1
1741  IF(j > 1) mj=j-india(j)+india(j-1)+1 ! first index in row J
1742  kj=india(k)-k+j ! index kj
1743  diag=c(india(j)) ! j-th diagonal element
1744 
1745  DO i=max(mj,mk),j-1
1746  ! L_kj = L_kj - L_ki *D_ii *L_ji
1747  c(kj)=c(kj) - c(india(k)-k+i)*c(india(j)-j+i)
1748  END DO ! I
1749 
1750  IF(j /= k) c(kj)=c(kj)*diag
1751  END DO ! J
1752 
1753  IF(diag+c(india(k)) > diag) THEN ! test for linear dependence
1754  c(india(k))=1.0d0/sqrt(c(india(k))) ! square root
1755  ELSE
1756  DO j=mk,k ! reset row K
1757  c(india(k)-k+j)=0.0d0
1758  END DO ! J
1759  IF(nrkd == 0) THEN
1760  nrkd=-k
1761  ELSE
1762  IF(nrkd < 0) nrkd=1
1763  nrkd=nrkd+1
1764  END IF
1765  END IF
1766 
1767  END DO ! K
1768  return
1769 END SUBROUTINE lltdec
1770 
1771 
1783 
1784 SUBROUTINE lltfwd(n,c,india,x)
1785  IMPLICIT NONE
1786  INTEGER :: j
1787  INTEGER :: k
1788 
1789  INTEGER, INTENT(IN) :: n
1790  DOUBLE PRECISION, INTENT(IN) :: c(*)
1791  INTEGER, INTENT(IN) :: india(n)
1792  DOUBLE PRECISION, INTENT(IN OUT) :: x(n)
1793 
1794  x(1)=x(1)*c(india(1))
1795  DO k=2,n ! forward loop
1796  DO j=k-india(k)+india(k-1)+1,k-1
1797  x(k)=x(k)-c(india(k)-k+j)*x(j) ! X_k := X_k - L_kj * B_j
1798  END DO ! J
1799  x(k)=x(k)*c(india(k))
1800  END DO ! K
1801  return
1802 END SUBROUTINE lltfwd
1803 
1815 
1816 SUBROUTINE lltbwd(n,c,india,x)
1817  IMPLICIT NONE
1818  INTEGER :: j
1819  INTEGER :: k
1820 
1821  INTEGER, INTENT(IN) :: n
1822  DOUBLE PRECISION, INTENT(IN) :: c(*)
1823  INTEGER, INTENT(IN) :: india(n)
1824  DOUBLE PRECISION, INTENT(IN OUT) :: x(n)
1825 
1826  DO k=n,2,-1 ! backward loop
1827  x(k)=x(k)*c(india(k))
1828  DO j=k-india(k)+india(k-1)+1,k-1
1829  x(j)=x(j)-c(india(k)-k+j)*x(k) ! X_j := X_j - L_kj * X_k
1830  END DO ! J
1831  END DO ! K
1832  x(1)=x(1)*c(india(1))
1833 END SUBROUTINE lltbwd
1834 
1851 SUBROUTINE equdec(n,m,c,india,nrkd,nrkd2)
1852  IMPLICIT NONE
1853  REAL:: eps
1854  INTEGER:: i
1855  INTEGER:: j
1856  INTEGER:: jk
1857  INTEGER:: k
1858  INTEGER:: ntotal
1859 
1860  INTEGER, INTENT(IN) :: n
1861  INTEGER, INTENT(IN) :: m
1862  DOUBLE PRECISION, INTENT(IN OUT) :: c(*)
1863  INTEGER, INTENT(IN OUT) :: india(n+m)
1864  INTEGER, INTENT(OUT) :: nrkd
1865  INTEGER, INTENT(OUT) :: nrkd2
1866 
1867 
1868  parameter(eps=2.22044605e-16)
1869  ! ...
1870  ntotal=n+n*m+(m*m+m)/2
1871 
1872  CALL lltdec(n,c,india,nrkd) ! decomposition G G^T
1873  DO i=1,m
1874  CALL lltfwd(n,c,india,c(india(n)+(i-1)*n+1)) ! forward solution K
1875  END DO
1876 
1877  jk=india(n)+n*m
1878  DO j=1,m
1879  DO k=1,j
1880  jk=jk+1
1881  c(jk)=0.0d0 ! product K K^T
1882  DO i=1,n
1883  c(jk)=c(jk)+c(india(n)+(j-1)*n+i)*c(india(n)+(k-1)*n+i)
1884  END DO
1885  END DO
1886  END DO
1887 
1888  india(n+1)=1
1889  DO i=2,m
1890  india(n+i)=india(n+i-1)+min(i,m) ! pointer for K K^T
1891  END DO
1892 
1893  CALL lltdec(m,c(india(n)+n*m+1),india(n+1),nrkd2) ! decomp. H H^T
1894 
1895  ntotal=n+n*m+(m*m+m)/2
1896 
1897  return
1898  END SUBROUTINE equdec
1899 
1915  SUBROUTINE equslv(n,m,c,india,x) ! solution vector
1916  IMPLICIT NONE
1917  INTEGER :: i
1918  INTEGER :: j
1919 
1920  INTEGER, INTENT(IN) :: n
1921  INTEGER, INTENT(IN) :: m
1922  DOUBLE PRECISION, INTENT(IN) :: c(*)
1923  INTEGER, INTENT(IN) :: india(n+m)
1924  DOUBLE PRECISION, INTENT(IN OUT) :: x(n+m)
1925 
1926  CALL lltfwd(n,c,india,x) ! result is u
1927  DO i=1,m
1928  DO j=1,n
1929  x(n+i)=x(n+i)-x(j)*c(india(n)+(i-1)*n+j) ! g - K u
1930  END DO
1931  END DO
1932  CALL lltfwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is v
1933 
1934 
1935  CALL lltbwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is -y
1936  DO i=1,m
1937  x(n+i)=-x(n+i) ! result is +y
1938  END DO
1939 
1940  DO i=1,n
1941  DO j=1,m
1942  x(i)=x(i)-x(n+j)*c(india(n)+(j-1)*n+i) ! u - K^T y
1943  END DO
1944  END DO
1945  CALL lltbwd(n,c,india,x) ! result is x
1946 END SUBROUTINE equslv
1947 
1976 
1977 SUBROUTINE precon(p,n,c,cu,a,s)
1978  IMPLICIT NONE
1979  INTEGER :: i
1980  INTEGER :: ii
1981  INTEGER :: j
1982  INTEGER :: jj
1983  INTEGER :: jk
1984  INTEGER :: k
1985  INTEGER :: kk
1986 
1987  INTEGER, INTENT(IN) :: p
1988  INTEGER, INTENT(IN) :: n
1989  DOUBLE PRECISION, INTENT(IN) :: c(n)
1990  DOUBLE PRECISION, INTENT(OUT) :: cu(n)
1991  DOUBLE PRECISION, INTENT(IN OUT) :: a(n,p)
1992  DOUBLE PRECISION, INTENT(OUT) :: s((p*p+p)/2)
1993 
1994  DOUBLE PRECISION :: div
1995  DOUBLE PRECISION :: ratio
1996 
1997  DO i=1,(p*p+p)/2
1998  s(i)=0.0d0
1999  END DO
2000  DO i=1,n
2001  jk=0
2002  div=c(i) ! copy
2003  IF (div > 0.0d0) THEN
2004  cu(i)=1.0d0/dsqrt(div)
2005  ELSE
2006  cu(i)=0.0d0
2007  END IF
2008  DO j=1,p
2009  a(i,j)=a(i,j)*cu(i) ! K = A C^{-1/2}
2010  DO k=1,j
2011  jk=jk+1
2012  s(jk)=s(jk)+a(i,j)*a(i,k) ! S = symmetric matrix K K^T
2013  END DO
2014  END DO
2015  END DO
2016 
2017  ii=0
2018  DO i=1,p ! S -> H D H^T (Cholesky)
2019  ii=ii+i
2020  IF(s(ii) /= 0.0d0) s(ii)=1.0d0/s(ii)
2021  jj=ii
2022  DO j=i+1,p
2023  ratio=s(i+jj)*s(ii)
2024  kk=jj
2025  DO k=j,p
2026  s(kk+j)=s(kk+j)-s(kk+i)*ratio
2027  kk=kk+k
2028  END DO ! K
2029  s(i+jj)=ratio
2030  jj=jj+j
2031  END DO ! J
2032  END DO ! I
2033  return
2034  END SUBROUTINE precon
2035 
2045 
2046  SUBROUTINE presol(p,n,cu,a,s,x,y) ! solution
2047  IMPLICIT NONE
2048  INTEGER :: i
2049  INTEGER :: j
2050  INTEGER :: jj
2051  INTEGER :: k
2052  INTEGER :: kk
2053 
2054  INTEGER, INTENT(IN) :: p
2055  INTEGER, INTENT(IN) :: n
2056 
2057  DOUBLE PRECISION, INTENT(IN) :: cu(n)
2058  DOUBLE PRECISION, INTENT(IN) :: a(n,p)
2059  DOUBLE PRECISION, INTENT(IN) :: s((p*p+p)/2)
2060  DOUBLE PRECISION, INTENT(OUT) :: x(n+p)
2061  DOUBLE PRECISION, INTENT(IN) :: y(n+p)
2062 
2063  DOUBLE PRECISION :: dsum
2064 
2065  DO i=1,n+p
2066  x(i)=y(i)
2067  END DO
2068  DO i=1,n
2069  x(i)=x(i)*cu(i) ! u =C^{-1/2} y
2070  DO j=1,p
2071  x(n+j)=x(n+j)-a(i,j)*x(i) ! d - K u
2072  END DO
2073  END DO
2074 
2075  jj=0
2076  DO j=1,p ! Cholesky solution for v
2077  dsum=x(n+j)
2078  DO k=1,j-1
2079  dsum=dsum-s(k+jj)*x(n+k) ! H v = d - K u
2080  END DO
2081  x(n+j)=dsum ! -> v
2082  jj=jj+j
2083  END DO
2084 
2085  DO j=p,1,-1 ! solution for lambda
2086  dsum=x(n+j)*s(jj)
2087  kk=jj
2088  DO k=j+1,p
2089  dsum=dsum+s(kk+j)*x(n+k) ! D H^T lambda = -v
2090  kk=kk+k
2091  END DO
2092  x(n+j)=-dsum ! -> lambda
2093  jj=jj-j
2094  END DO
2095 
2096  DO i=1,n ! u - K^T lambda
2097  DO j=1,p
2098  x(i)=x(i)-a(i,j)*x(n+j)
2099  END DO
2100  END DO
2101  DO i=1,n
2102  x(i)=x(i)*cu(i) ! x = C^{-1/2} u
2103  END DO
2104 
2105 END SUBROUTINE presol
2106 
2107 
2108 ! 090817 C. Kleinwort, DESY-FH1
2153 SUBROUTINE sqmibb(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
2154  ! Double precision scratch arrays:
2155  ! VBND(N*(NBND+1)) = storage of band part
2156  ! VBDR(N* NBDR) = storage of border part
2157  ! AUX (N* NBDR) = intermediate results
2158 
2159  ! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only)
2160 
2161  IMPLICIT NONE
2162  INTEGER:: i
2163  INTEGER:: ib
2164  INTEGER:: ij
2165  INTEGER:: ioff
2166  INTEGER:: ip
2167  INTEGER:: ip1
2168  INTEGER:: ip2
2169  INTEGER:: is
2170  INTEGER:: j
2171  INTEGER:: j0
2172  INTEGER:: jb
2173  INTEGER:: joff
2174  INTEGER:: mp1
2175  INTEGER:: nb1
2176  INTEGER:: nmb
2177  INTEGER:: npri
2178  INTEGER:: nrankb
2179 
2180  DOUBLE PRECISION, INTENT(IN OUT) :: v(*)
2181  DOUBLE PRECISION, INTENT(OUT) :: b(n)
2182  INTEGER, INTENT(IN) :: n
2183  INTEGER, INTENT(IN) :: nbdr
2184  INTEGER, INTENT(IN) :: nbnd
2185  INTEGER, INTENT(IN) :: inv
2186  INTEGER, INTENT(OUT) :: nrank
2187 
2188  DOUBLE PRECISION, INTENT(OUT) :: vbnd(n*(nbnd+1))
2189  DOUBLE PRECISION, INTENT(OUT) :: vbdr(n*nbdr)
2190  DOUBLE PRECISION, INTENT(OUT) :: aux(n*nbdr)
2191  DOUBLE PRECISION, INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
2192  DOUBLE PRECISION, INTENT(OUT) :: vzru(nbdr)
2193  DOUBLE PRECISION, INTENT(OUT) :: scdiag(nbdr)
2194  INTEGER, INTENT(OUT) :: scflag(nbdr)
2195 
2196  SAVE npri
2197  DATA npri / 100 /
2198  ! ...
2199  nrank=0
2200  nb1=nbdr+1
2201  mp1=nbnd+1
2202  nmb=n-nbdr
2203  ! copy band part
2204  DO i=nb1,n
2205  ip=(i*(i+1))/2
2206  is=0
2207  DO j=i,min(n,i+nbnd)
2208  ip=ip+is
2209  is=j
2210  ib=j-i+1
2211  vbnd(ib+(i-nb1)*mp1)=v(ip)
2212  END DO
2213  END DO
2214  ! copy border part
2215  IF (nbdr > 0) THEN
2216  ioff=0
2217  DO i=1,nbdr
2218  ip=(i*(i+1))/2
2219  is=0
2220  DO j=i,n
2221  ip=ip+is
2222  is=j
2223  vbdr(ioff+j)=v(ip)
2224  END DO
2225  ioff=ioff+n
2226  END DO
2227  END IF
2228 
2229  CALL dbcdec(vbnd,mp1,nmb,aux)
2230  ! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable
2231  ! CALL DBCPRB(VBND,MP1,NMB)
2232  ip=1
2233  DO i=1, nmb
2234  IF (vbnd(ip) <= 0.0d0) THEN
2235  npri=npri-1
2236  IF (npri >= 0) THEN
2237  IF (vbnd(ip) == 0.0d0) THEN
2238  print *, ' SQMIBB matrix singular', n, nbdr, nbnd
2239  ELSE
2240  print *, ' SQMIBB matrix not positive definite', n, nbdr, nbnd
2241  END IF
2242  END IF
2243  ! return zeros
2244  DO ip=1,n
2245  b(ip)=0.0d0
2246  END DO
2247  DO ip=1,(n*n+n)/2
2248  v(ip)=0.0d0
2249  END DO
2250  return
2251  END IF
2252  ip=ip+mp1
2253  END DO
2254  nrank=nmb
2255 
2256  IF (nbdr == 0) THEN ! special case NBDR=0
2257 
2258  CALL dbcslv(vbnd,mp1,nmb,b,b)
2259  IF (inv > 0) THEN
2260  IF (inv > 1) THEN
2261  CALL dbcinv(vbnd,mp1,nmb,v)
2262  ELSE
2263  CALL dbcinb(vbnd,mp1,nmb,v)
2264  END IF
2265  END IF
2266 
2267  ELSE ! general case NBDR>0
2268 
2269  ioff=nb1
2270  DO ib=1,nbdr
2271  ! solve for aux. vectors
2272  CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff),aux(ioff))
2273  ! zT ru
2274  vzru(ib)=b(ib)
2275  DO i=0,nmb-1
2276  vzru(ib)=vzru(ib)-b(nb1+i)*aux(ioff+i)
2277  END DO
2278  ioff=ioff+n
2279  END DO
2280  ! solve for band part only
2281  CALL dbcslv(vbnd,mp1,nmb,b(nb1),b(nb1))
2282  ! Ck - cT z
2283  ip=0
2284  ioff=nb1
2285  DO ib=1,nbdr
2286  joff=nb1
2287  DO jb=1,ib
2288  ip=ip+1
2289  vbk(ip)=v(ip)
2290  DO i=0,nmb-1
2291  vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
2292  END DO
2293  joff=joff+n
2294  END DO
2295  ioff=ioff+n
2296  END DO
2297  ! solve border part
2298  CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
2299  IF (nrankb == nbdr) THEN
2300  nrank=nrank+nbdr
2301  ELSE
2302  npri=npri-1
2303  IF (npri >= 0) print *, ' SQMIBB undef border ', n, nbdr, nbnd, nrankb
2304  DO ib=1,nbdr
2305  vzru(ib)=0.0d0
2306  END DO
2307  DO ip=(nbdr*nbdr+nbdr)/2,1,-1
2308  vbk(ip)=0.0d0
2309  END DO
2310  END IF
2311  ! smoothed data points
2312  ioff=nb1
2313  DO ib=1, nbdr
2314  b(ib) = vzru(ib)
2315  DO i=0,nmb-1
2316  b(nb1+i)=b(nb1+i)-b(ib)*aux(ioff+i)
2317  END DO
2318  ioff=ioff+n
2319  END DO
2320  ! inverse requested ?
2321  IF (inv > 0) THEN
2322  IF (inv > 1) THEN
2323  CALL dbcinv(vbnd,mp1,nmb,v)
2324  ELSE
2325  CALL dbcinb(vbnd,mp1,nmb,v)
2326  END IF
2327  ! expand/correct from NMB to N
2328  ip1=(nmb*nmb+nmb)/2
2329  ip2=(n*n+n)/2
2330  DO i=nmb-1,0,-1
2331  j0=0
2332  IF (inv == 1) j0=max(0,i-nbnd)
2333  DO j=i,j0,-1
2334  v(ip2)=v(ip1)
2335  ioff=nb1
2336  DO ib=1,nbdr
2337  joff=nb1
2338  DO jb=1,nbdr
2339  ij=max(ib,jb)
2340  ij=(ij*ij-ij)/2+min(ib,jb)
2341  v(ip2)=v(ip2)+vbk(ij)*aux(ioff+i)*aux(joff+j)
2342  joff=joff+n
2343  END DO
2344  ioff=ioff+n
2345  END DO
2346  ip1=ip1-1
2347  ip2=ip2-1
2348  END DO
2349  ip1=ip1-j0
2350  ip2=ip2-j0
2351 
2352  DO ib=nbdr,1,-1
2353  v(ip2)=0.0d0
2354  joff=nb1
2355  DO jb=1,nbdr
2356  ij=max(ib,jb)
2357  ij=(ij*ij-ij)/2+min(ib,jb)
2358  v(ip2)=v(ip2)-vbk(ij)*aux(i+joff)
2359  joff=joff+n
2360  END DO
2361  ip2=ip2-1
2362  END DO
2363  END DO
2364 
2365  DO ip=(nbdr*nbdr+nbdr)/2,1,-1
2366  v(ip2)=vbk(ip)
2367  ip2=ip2-1
2368  END DO
2369 
2370  END IF
2371  END IF
2372 
2373 END SUBROUTINE sqmibb
Generated on Wed May 23 2012 16:53:05 for Millepede-II by   doxygen 1.8.1