Millepede-II  V04-08-03
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 
75 
76 !----------------------------------------------------------------------
95 
96 SUBROUTINE sqminv(v,b,n,nrank,diag,next) ! matrix inversion
97  USE mpdef
98 
99  IMPLICIT NONE
100  INTEGER(mpi) :: i
101  INTEGER(mpi) :: ij
102  INTEGER(mpi) :: j
103  INTEGER(mpi) :: jj
104  INTEGER(mpi) :: jk
105  INTEGER(mpi) :: jl
106  INTEGER(mpi) :: k
107  INTEGER(mpi) :: kk
108  INTEGER(mpi) :: l
109  INTEGER(mpi) :: last
110  INTEGER(mpi) :: lk
111  INTEGER(mpi) :: next0
112 
113  REAL(mpd), INTENT(IN OUT) :: v(*)
114  REAL(mpd), INTENT(OUT) :: b(n)
115  INTEGER(mpi), INTENT(IN) :: n
116  INTEGER(mpi), INTENT(OUT) :: nrank
117  REAL(mpd), INTENT(OUT) :: diag(n)
118  INTEGER(mpi), INTENT(OUT) :: next(n)
119  REAL(mpd) :: vkk
120  REAL(mpd) :: vjk
121 
122  !REAL(mpd), PARAMETER :: eps=1.0E-10_mpd
123  REAL(mpd) eps
124  ! ...
125  eps = 16.0_mpd * epsilon(eps) ! 16 * precision(mpd)
126 
127  next0=1
128  l=1
129  DO i=1,n
130  next(i)=i+1 ! set "next" pointer
131  diag(i)=abs(v((i*i+i)/2)) ! save abs of diagonal elements
132  END DO
133  next(n)=-1 ! end flag
134 
135  nrank=0
136  DO i=1,n ! start of loop
137  k =0
138  vkk=0.0_mpd
139 
140  j=next0
141  last=0
142 05 IF(j > 0) THEN
143  jj=(j*j+j)/2
144  IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
145  vkk=v(jj)
146  k=j
147  l=last
148  END IF
149  last=j
150  j=next(last)
151  GO TO 05
152  END IF
153 
154  IF(k /= 0) THEN ! pivot found
155  kk=(k*k+k)/2
156  IF(l == 0) THEN
157  next0=next(k)
158  ELSE
159  next(l)=next(k)
160  END IF
161  next(k)=0 ! index is used, reset
162  nrank=nrank+1 ! increase rank and ...
163  vkk =1.0/vkk
164  v(kk) =-vkk
165  b(k) =b(k)*vkk
166  jk =kk-k
167  jl =0
168  DO j=1,n ! elimination
169  IF(j == k) THEN
170  jk=kk
171  jl=jl+j
172  ELSE
173  IF(j < k) THEN
174  jk=jk+1
175  ELSE
176  jk=jk+j-1
177  END IF
178  vjk =v(jk)
179  v(jk)=vkk*vjk
180  b(j) =b(j)-b(k)*vjk
181  lk =kk-k
182  DO l=1,j
183  jl=jl+1
184  IF(l == k) THEN
185  lk=kk
186  ELSE
187  IF(l < k) THEN
188  lk=lk+1
189  ELSE
190  lk=lk+l-1
191  END IF
192  v(jl)=v(jl)-v(lk)*vjk
193  END IF
194  END DO
195  END IF
196  END DO
197  ELSE
198  DO k=1,n
199  IF(next(k) /= 0) THEN
200  b(k)=0.0_mpd ! clear vector element
201  DO j=1,k
202  IF(next(j) /= 0) v((k*k-k)/2+j)=0.0_mpd ! clear matrix row/col
203  END DO
204  END IF
205  END DO
206  GO TO 10
207  END IF
208  END DO ! end of loop
209  10 DO ij=1,(n*n+n)/2
210  v(ij)=-v(ij) ! finally reverse sign of all matrix elements
211  END DO
212 END SUBROUTINE sqminv
213 
228 
229 SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk,mon) !
230  USE mpdef
231 
232  IMPLICIT NONE
233  INTEGER(mpi) :: i
234  INTEGER(mpi) :: j
235  INTEGER(mpi) :: k
236  INTEGER(mpi) :: l
237  INTEGER(mpi) :: last
238  INTEGER(mpi) :: next0
239 
240  REAL(mpd), INTENT(IN OUT) :: v(*)
241  REAL(mpd), INTENT(OUT) :: b(n)
242  INTEGER(mpi), INTENT(IN) :: n
243  INTEGER(mpi), INTENT(OUT) :: nrank
244  REAL(mpd), INTENT(OUT) :: diag(n)
245  INTEGER(mpi), INTENT(OUT) :: next(n)
246  REAL(mpd), INTENT(OUT) :: vk(n)
247  INTEGER(mpi), INTENT(IN) :: mon
248 
249  INTEGER(mpl) :: i8
250  INTEGER(mpl) :: j8
251  INTEGER(mpl) :: jj
252  INTEGER(mpl) :: k8
253  INTEGER(mpl) :: kk
254  INTEGER(mpl) :: kkmk
255  INTEGER(mpl) :: jk
256  INTEGER(mpl) :: jl
257 
258  REAL(mpd) :: vkk
259  REAL(mpd) :: vjk
260 
261  REAL(mpd), PARAMETER :: eps=1.0e-10_mpd
262  ! ...
263  next0=1
264  l=1
265  DO i=1,n
266  i8=int(i,mpl)
267  next(i)=i+1 ! set "next" pointer
268  diag(i)=abs(v((i8*i8+i8)/2)) ! save abs of diagonal elements
269  END DO
270  next(n)=-1 ! end flag
271 
272  nrank=0
273  DO i=1,n ! start of loop
274  ! monitoring ?
275  IF(mon>0) CALL monpgs(i)
276  k =0
277  vkk=0.0_mpd
278  j=next0
279  last=0
280 05 IF(j > 0) THEN
281  j8=int(j,mpl)
282  jj=(j8*j8+j8)/2
283  IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
284  vkk=v(jj)
285  k=j
286  l=last
287  END IF
288  last=j
289  j=next(last)
290  GO TO 05
291  END IF
292 
293  IF(k /= 0) THEN ! pivot found
294  k8=int(k,mpl)
295  kk=(k8*k8+k8)/2
296  kkmk=kk-k8
297  IF(l == 0) THEN
298  next0=next(k)
299  ELSE
300  next(l)=next(k)
301  END IF
302  next(k)=0 ! index is used, reset
303  nrank=nrank+1 ! increase rank and ...
304  vkk =1.0/vkk
305  v(kk) =-vkk
306  b(k) =b(k)*vkk
307  ! elimination
308  jk =kkmk
309  DO j=1,n
310  IF(j == k) THEN
311  jk=kk
312  vk(j)=0.
313  ELSE
314  IF(j < k) THEN
315  jk=jk+1
316  ELSE
317  jk=jk+int(j,mpl)-1
318  END IF
319  v(jk)=v(jk)*vkk
320  vk(j)=v(jk)
321  END IF
322  END DO
323  ! parallelize row loop
324  ! slot of 128 'J' for next idle thread (optimized on Intel Xeon)
325  !$OMP PARALLEL DO &
326  !$OMP PRIVATE(JL,VJK,J8) &
327  !$OMP SCHEDULE(DYNAMIC,128)
328  DO j=n,1,-1
329  IF(j == k) cycle
330  j8=int(j,mpl)
331  jl=j8*(j8-1)/2
332  vjk =vk(j)/vkk
333  b(j) =b(j)-b(k)*vjk
334  v(jl+1:jl+j)=v(jl+1:jl+j)-vk(1:j)*vjk
335  END DO
336  !$OMP END PARALLEL DO
337  ELSE
338  DO k=1,n
339  k8=int(k,mpl)
340  kk=(k8*k8-k8)/2
341  IF(next(k) /= 0) THEN
342  b(k)=0.0_mpd ! clear vector element
343  DO j=1,k
344  IF(next(j) /= 0) v(kk+int(j,mpl))=0.0_mpd ! clear matrix row/col
345  END DO
346  END IF
347  END DO
348  GO TO 10
349  END IF
350  END DO ! end of loop
351  10 DO jj=1,(int(n,mpl)*int(n+1,mpl))/2
352  v(jj)=-v(jj) ! finally reverse sign of all matrix elements
353  END DO
354 END SUBROUTINE sqminl
355 
367 
368 SUBROUTINE devrot(n,diag,u,v,work,iwork) ! diagonalization
369  USE mpdef
370 
371  IMPLICIT NONE
372 
373  INTEGER(mpi), INTENT(IN) :: n
374  REAL(mpd), INTENT(OUT) :: diag(n)
375  REAL(mpd), INTENT(OUT) :: u(n,n)
376  REAL(mpd), INTENT(IN) :: v(*)
377  REAL(mpd), INTENT(OUT) :: work(n)
378  INTEGER(mpi), INTENT(OUT) :: iwork(n)
379 
380 
381  INTEGER(mpi), PARAMETER :: itmax=30
382  REAL(mpd), PARAMETER :: tol=epsilon(tol)
383  REAL(mpd), PARAMETER :: eps=epsilon(eps)
384 
385  REAL(mpd) :: f
386  REAL(mpd) :: g
387  REAL(mpd) :: h
388  REAL(mpd) :: sh
389  REAL(mpd) :: hh
390  REAL(mpd) :: b
391  REAL(mpd) :: p
392  REAL(mpd) :: r
393  REAL(mpd) :: s
394  REAL(mpd) :: c
395  REAL(mpd) :: workd
396 
397  INTEGER(mpi) :: ij
398  INTEGER(mpi) :: i
399  INTEGER(mpi) :: j
400  INTEGER(mpi) :: k
401  INTEGER(mpi) :: l
402  INTEGER(mpi) :: m
403  INTEGER(mpi) :: ll
404  ! ...
405  ! 1. part: symmetric matrix V reduced to tridiagonal from
406  ij=0
407  DO i=1,n
408  DO j=1,i
409  ij=ij+1
410  u(i,j)=v(ij) ! copy half of symmetric matirx
411  END DO
412  END DO
413 
414  DO i=n,2,-1
415  l=i-2
416  f=u(i,i-1)
417  g=0.0_mpd
418  IF(l /= 0) THEN
419  DO k=1,l
420  IF(abs(u(i,k)) > tol) g=g+u(i,k)*u(i,k)
421  END DO
422  h=g+f*f
423  END IF
424  IF(g < tol) THEN ! G too small
425  work(i)=f ! skip transformation
426  h =0.0_mpd
427  ELSE
428  l=l+1
429  sh=sqrt(h)
430  IF(f >= 0.0_mpd) sh=-sh
431  g=sh
432  work(i)=sh
433  h=h-f*g
434  u(i,i-1)=f-g
435  f=0.0_mpd
436  DO j=1,l
437  u(j,i)=u(i,j)/h
438  g=0.0_mpd
439  ! form element of a u
440  DO k=1,j
441  IF(abs(u(j,k)) > tol.AND.abs(u(i,k)) > tol) THEN
442  g=g+u(j,k)*u(i,k)
443  END IF
444  END DO
445  DO k=j+1,l
446  IF(abs(u(k,j)) > tol.AND.abs(u(i,k)) > tol) THEN
447  g=g+u(k,j)*u(i,k)
448  END IF
449  END DO
450  work(j)=g/h
451  f=f+g*u(j,i)
452  END DO
453  ! form k
454  hh=f/(h+h)
455  ! form reduced a
456  DO j=1,l
457  f=u(i,j)
458  work(j)=work(j)-hh*f
459  g=work(j)
460  DO k=1,j
461  u(j,k)=u(j,k)-f*work(k)-g*u(i,k)
462  END DO
463  END DO
464  END IF
465  diag(i)=h
466  END DO
467 
468  diag(1)=0.0_mpd
469  work(1)=0.0_mpd
470 
471  ! accumulation of transformation matrices
472  DO i=1,n
473  IF(diag(i) /= 0.0) THEN
474  DO j=1,i-1
475  g=0.0_mpd
476  DO k=1,i-1
477  g=g+u(i,k)*u(k,j)
478  END DO
479  DO k=1,i-1
480  u(k,j)=u(k,j)-g*u(k,i)
481  END DO
482  END DO
483  END IF
484  diag(i)=u(i,i)
485  u(i,i)=1.0_mpd
486  DO j=1,i-1
487  u(i,j)=0.0_mpd
488  u(j,i)=0.0_mpd
489  END DO
490  END DO
491 
492  ! 2. part: diagonalization of tridiagonal matrix
493  DO i=2,n
494  work(i-1)=work(i)
495  END DO
496  work(n)=0.0_mpd
497  b=0.0_mpd
498  f=0.0_mpd
499 
500  DO l=1,n
501  j=0
502  h=eps*(abs(diag(l))+abs(work(l)))
503  IF(b < h) b=h
504  DO m=l,n
505  IF(abs(work(m)) <= b) GO TO 10 ! look for small sub-diagonal element
506  END DO
507  m=l
508 10 IF(m == l) GO TO 30
509  ! next iteration
510 20 IF(j == itmax) THEN
511  WRITE(*,*) 'DEVROT: Iteration limit reached'
512  CALL peend(32,'Aborted, iteration limit reached in diagonalization')
513  stop
514  END IF
515  j=j+1
516  g=diag(l)
517  p=(diag(l+1)-g)/(2.0_mpd*work(l))
518  r=sqrt(1.0_mpd+p*p)
519  diag(l)=work(l)
520  IF(p < 0.0_mpd) diag(l)=diag(l)/(p-r)
521  IF(p >= 0.0_mpd) diag(l)=diag(l)/(p+r)
522  h=g-diag(l)
523  DO i=l+1,n
524  diag(i)=diag(i)-h
525  END DO
526  f=f+h
527  ! QL transformation
528  p=diag(m)
529  c=1.0_mpd
530  s=0.0_mpd
531  DO i=m-1,l,-1 ! reverse loop
532  g=c*work(i)
533  h=c*p
534  IF(abs(p) >= abs(work(i))) THEN
535  c=work(i)/p
536  r=sqrt(1.0_mpd+c*c)
537  work(i+1)=s*p*r
538  s=c/r
539  c=1.0_mpd/r
540  ELSE
541  c=p/work(i)
542  r=sqrt(1.0_mpd+c*c)
543  work(i+1)=s*work(i)*r
544  s=1.0_mpd/r
545  c=c/r
546  END IF
547  p=c*diag(i)-s*g
548  diag(i+1)=h+s*(c*g+s*diag(i))
549  ! form vector
550  DO k=1,n
551  h=u(k,i+1)
552  u(k,i+1)=s*u(k,i)+c*h
553  u(k,i)=c*u(k,i)-s*h
554  END DO
555  END DO
556  work(l)=s*p
557  diag(l)=c*p
558  IF(abs(work(l)) > b) GO TO 20 ! next iteration
559 30 diag(l)=diag(l)+f
560  END DO
561  DO i=1,n
562  iwork(i)=i
563  END DO
564 
565  m=1
566 40 m=1+3*m ! determine initial increment
567  IF(m <= n) GO TO 40
568 50 m=m/3
569  DO j=1,n-m ! sort with increment M
570  l=j
571 60 IF(diag(iwork(l+m)) > diag(iwork(l))) THEN ! compare
572  ll=iwork(l+m) ! exchange the two index values
573  iwork(l+m)=iwork(l)
574  iwork(l)=ll
575  l=l-m
576  IF(l > 0) GO TO 60
577  END IF
578  END DO
579  IF(m > 1) GO TO 50
580 
581  DO i=1,n
582  IF(iwork(i) /= i) THEN
583  ! move vector from position I to the work area
584  workd=diag(i)
585  DO l=1,n
586  work(l)=u(l,i)
587  END DO
588  k=i
589 70 j=k
590  k=iwork(j)
591  iwork(j)=j
592  IF(k /= i) THEN
593  ! move vector from position K to the (free) position J
594  diag(j)=diag(k)
595  DO l=1,n
596  u(l,j)=u(l,k)
597  END DO
598  GO TO 70
599  END IF
600  ! move vector from the work area to position J
601  diag(j)=workd
602  DO l=1,n
603  u(l,j)=work(l)
604  END DO
605  END IF
606  END DO
607 END SUBROUTINE devrot
608 
610 SUBROUTINE devsig(n,diag,u,b,coef)
611  USE mpdef
612 
613  IMPLICIT NONE
614 
615  INTEGER(mpi), INTENT(IN) :: n
616  REAL(mpd), INTENT(IN) :: diag(n)
617  REAL(mpd), INTENT(IN) :: u(n,n)
618  REAL(mpd), INTENT(IN) :: b(n)
619  REAL(mpd), INTENT(OUT) :: coef(n)
620  INTEGER(mpi) :: i
621  INTEGER(mpi) :: j
622  REAL(mpd) :: dsum
623  ! ...
624  DO i=1,n
625  coef(i)=0.0_mpd
626  IF(diag(i) > 0.0_mpd) THEN
627  dsum=0.0_mpd
628  DO j=1,n
629  dsum=dsum+u(j,i)*b(j)
630  END DO
631  coef(i)=abs(dsum)/sqrt(diag(i))
632  END IF
633  END DO
634 END SUBROUTINE devsig
635 
636 
647 
648 SUBROUTINE devsol(n,diag,u,b,x,work)
649  USE mpdef
650 
651  IMPLICIT NONE
652 
653  INTEGER(mpi), INTENT(IN) :: n
654  REAL(mpd), INTENT(IN) :: diag(n)
655  REAL(mpd), INTENT(IN) :: u(n,n)
656  REAL(mpd), INTENT(IN) :: b(n)
657  REAL(mpd), INTENT(OUT) :: x(n)
658  REAL(mpd), INTENT(OUT) :: work(n)
659  INTEGER(mpi) :: i
660  INTEGER(mpi) :: j
661  INTEGER(mpi) :: jj
662  REAL(mpd) :: s
663  ! ...
664  DO j=1,n
665  s=0.0_mpd
666  work(j)=0.0_mpd
667  IF(diag(j) /= 0.0_mpd) THEN
668  DO i=1,n
669  ! j-th eigenvector is U(.,J)
670  s=s+u(i,j)*b(i)
671  END DO
672  work(j)=s/diag(j)
673  END IF
674  END DO
675 
676  DO j=1,n
677  s=0.0_mpd
678  DO jj=1,n
679  s=s+u(j,jj)*work(jj)
680  END DO
681  x(j)=s
682  END DO
683 ! WRITE(*,*) 'DEVSOL'
684 ! WRITE(*,*) 'X ',X
685 END SUBROUTINE devsol
686 
694 
695 SUBROUTINE devinv(n,diag,u,v)
696  USE mpdef
697 
698  IMPLICIT NONE
699  INTEGER(mpi) :: i
700  INTEGER(mpi) :: ij
701  INTEGER(mpi) :: j
702  INTEGER(mpi) :: k
703 
704  INTEGER(mpi), INTENT(IN) :: n
705  REAL(mpd), INTENT(IN) :: diag(n)
706  REAL(mpd), INTENT(IN) :: u(n,n)
707  REAL(mpd), INTENT(OUT) :: v(*)
708  REAL(mpd) :: dsum
709  ! ...
710  ij=0
711  DO i=1,n
712  DO j=1,i
713  ij=ij+1
714  dsum=0.0_mpd
715  DO k=1,n
716  IF(diag(k) /= 0.0_mpd) THEN
717  dsum=dsum+u(i,k)*u(j,k)/diag(k)
718  END IF
719  END DO
720  v(ij)=dsum
721  END DO
722  END DO
723 END SUBROUTINE devinv
724 
725 
744 SUBROUTINE choldc(g,n)
745  USE mpdef
746 
747  IMPLICIT NONE
748  INTEGER(mpi) :: i
749  INTEGER(mpi) :: ii
750  INTEGER(mpi) :: j
751  INTEGER(mpi) :: jj
752  INTEGER(mpi) :: k
753  INTEGER(mpi) :: kk
754 
755  REAL(mpd), INTENT(IN OUT) :: g(*)
756  INTEGER(mpi), INTENT(IN) :: n
757 
758  REAL(mpd) :: ratio
759  ! ...
760  ii=0
761  DO i=1,n
762  ii=ii+i
763  IF(g(ii) /= 0.0) g(ii)=1.0/g(ii) ! (I,I) div !
764  jj=ii
765  DO j=i+1,n
766  ratio=g(i+jj)*g(ii) ! (I,J) (I,I)
767  kk=jj
768  DO k=j,n
769  g(kk+j)=g(kk+j)-g(kk+i)*ratio ! (K,J) (K,I)
770  kk=kk+k
771  END DO ! K
772  g(i+jj)=ratio ! (I,J)
773  jj=jj+j
774  END DO ! J
775  END DO ! I
776  RETURN
777 END SUBROUTINE choldc
778 
790 SUBROUTINE cholsl(g,x,n)
791  USE mpdef
792 
793  IMPLICIT NONE
794  REAL(mpd) :: dsum
795  INTEGER(mpi) :: i
796  INTEGER(mpi) :: ii
797  INTEGER(mpi) :: k
798  INTEGER(mpi) :: kk
799 
800  REAL(mpd), INTENT(IN) :: g(*)
801  REAL(mpd), INTENT(IN OUT) :: x(n)
802  INTEGER(mpi), INTENT(IN) :: n
803 
804  ii=0
805  DO i=1,n
806  dsum=x(i)
807  DO k=1,i-1
808  dsum=dsum-g(k+ii)*x(k) ! (K,I)
809  END DO
810  x(i)=dsum
811  ii=ii+i
812  END DO
813  DO i=n,1,-1
814  dsum=x(i)*g(ii) ! (I,I)
815  kk=ii
816  DO k=i+1,n
817  dsum=dsum-g(kk+i)*x(k) ! (K,I)
818  kk=kk+k
819  END DO
820  x(i)=dsum
821  ii=ii-i
822  END DO
823  RETURN
824 END SUBROUTINE cholsl
825 
836 SUBROUTINE cholin(g,v,n)
837  USE mpdef
838 
839  IMPLICIT NONE
840  REAL(mpd) :: dsum
841  INTEGER(mpi) :: i
842  INTEGER(mpi) :: ii
843  INTEGER(mpi) :: j
844  INTEGER(mpi) :: k
845  INTEGER(mpi) :: l
846  INTEGER(mpi) :: m
847 
848  REAL(mpd), INTENT(IN) :: g(*)
849  REAL(mpd), INTENT( OUT) :: v(*)
850  INTEGER(mpi), INTENT(IN) :: n
851 
852  ii=(n*n-n)/2
853  DO i=n,1,-1
854  dsum=g(ii+i) ! (I,I)
855  DO j=i,1,-1
856  DO k=j+1,n
857  l=min(i,k)
858  m=max(i,k)
859  dsum=dsum-g(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K)
860  END DO
861  v(ii+j)=dsum ! (I,J)
862  dsum=0.0_mpd
863  END DO
864  ii=ii-i+1
865  END DO
866 END SUBROUTINE cholin
867 
868 ! 201026 C. Kleinwort, DESY-BELLE
890 SUBROUTINE chdec2(g,n,nrank,evmax,evmin,mon)
891  USE mpdef
892 
893  IMPLICIT NONE
894  INTEGER(mpi) :: i
895  INTEGER(mpi) :: j
896 
897  INTEGER(mpl) :: ii
898  INTEGER(mpl) :: jj
899  REAL(mpd) :: ratio
900 
901  REAL(mpd), INTENT(IN OUT) :: g(*)
902  INTEGER(mpi), INTENT(IN) :: n
903  INTEGER(mpi), INTENT(OUT) :: nrank
904  REAL(mpd), INTENT(OUT) :: evmin
905  REAL(mpd), INTENT(OUT) :: evmax
906  INTEGER(mpi), INTENT(IN) :: mon
907 
908  nrank=0
909  ii=(int(n,mpl)*int(n+1,mpl))/2
910  DO i=n,1,-1
911  ! monitoring ?
912  IF(mon>0) CALL monpgs(n+1-i)
913  IF (g(ii) > 0.0_mpd) THEN
914  ! update rank, min, max eigenvalue
915  nrank=nrank+1
916  IF (nrank == 1) THEN
917  evmax=g(ii)
918  evmin=g(ii)
919  ELSE
920  evmax=max(evmax,g(ii))
921  evmin=min(evmin,g(ii))
922  END IF
923  g(ii)=1.0/g(ii) ! (I,I) div !
924  END IF
925  ii=ii-i
926  ! parallelize row loop
927  ! slot of 32 'J' for next idle thread (optimized on Intel Xeon)
928  !$OMP PARALLEL DO &
929  !$OMP PRIVATE(RATIO,JJ) &
930  !$OMP SCHEDULE(DYNAMIC,32)
931  DO j=1,i-1
932  ratio=g(ii+j)*g(ii+i) ! (I,J) (I,I)
933  IF (ratio == 0.0_mpd) cycle
934  jj=(int(j-1,mpl)*int(j,mpl))/2
935  g(jj+1:jj+j)=g(jj+1:jj+j)-g(ii+1:ii+j)*ratio ! (K,J) (K,I)
936  END DO ! J
937  !$OMP END PARALLEL DO
938  g(ii+1:ii+i-1)=g(ii+1:ii+i-1)*g(ii+i) ! (I,J)
939  END DO ! I
940 
941 END SUBROUTINE chdec2
942 
943 ! 201026 C. Kleinwort, DESY-BELLE
952 SUBROUTINE chslv2(g,x,n)
953  USE mpdef
954 
955  IMPLICIT NONE
956  INTEGER(mpi) :: i
957  INTEGER(mpi) :: k
958  INTEGER(mpl) :: ii
959  INTEGER(mpl) :: kk
960  REAL(mpd) :: dsum
961 
962  REAL(mpd), INTENT(IN) :: g(*)
963  REAL(mpd), INTENT(IN OUT) :: x(n)
964  INTEGER(mpi), INTENT(IN) :: n
965 
966  ii=(int(n,mpl)*int(n+1,mpl))/2
967  DO i=n,1,-1
968  dsum=x(i)
969  kk=ii
970  DO k=i+1,n
971  dsum=dsum-g(kk+i)*x(k) ! (K,I)
972  kk=kk+k
973  END DO
974  x(i)=dsum
975  ii=ii-i
976  END DO
977  DO i=1,n
978  dsum=x(i)*g(ii+i) ! (I,I)
979  DO k=1,i-1
980  dsum=dsum-g(k+ii)*x(k) ! (K,I)
981  END DO
982  x(i)=dsum
983  ii=ii+i
984  END DO
985 
986 END SUBROUTINE chslv2
987 
988 ! variable band matrix operations ----------------------------------
989 
1010 
1011 SUBROUTINE vabdec(n,val,ilptr)
1012  USE mpdef
1013 
1014  IMPLICIT NONE
1015 
1016  INTEGER(mpi) :: i
1017  INTEGER(mpi) :: in
1018  INTEGER(mpi) :: j
1019  INTEGER(mpi) :: k
1020  INTEGER(mpi) :: kj
1021  INTEGER(mpi) :: mj
1022  INTEGER(mpi) :: mk
1023  REAL(mpd) :: sn
1024  REAL(mpd) :: beta
1025  REAL(mpd) :: delta
1026  REAL(mpd) :: theta
1027 
1028  INTEGER(mpi), INTENT(IN) :: n
1029  REAL(mpd), INTENT(IN OUT) :: val(*)
1030  INTEGER(mpi), INTENT(IN) :: ilptr(n)
1031 
1032  REAL(mpd) :: dgamma
1033  REAL(mpd) :: xi
1034  REAL(mpd) :: valkj
1035 
1036  REAL(mpd), PARAMETER :: one=1.0_mpd
1037  REAL(mpd), PARAMETER :: two=2.0_mpd
1038  REAL(mpd), PARAMETER :: eps = epsilon(eps)
1039 
1040  WRITE(*,*) 'Variable band matrix Cholesky decomposition'
1041 
1042  dgamma=0.0_mpd
1043  i=1
1044  DO j=1,ilptr(n) ! loop thrugh all matrix elements
1045  IF(ilptr(i) == j) THEN ! diagonal element
1046  IF(val(j) <= 0.0_mpd) GO TO 01 ! exit loop for negative diag
1047  dgamma=max(dgamma,abs(val(j))) ! max diagonal element
1048  i=i+1
1049  END IF
1050  END DO
1051  i=n+1
1052 01 in=i-1 ! IN positive diagonal elements
1053  WRITE(*,*) ' ',in,' positive diagonal elements'
1054  xi=0.0_mpd
1055  i=1
1056  DO j=1,ilptr(in) ! loop for positive diagonal elements
1057  ! through all matrix elements
1058  IF(ilptr(i) == j) THEN ! diagonal element
1059  i=i+1
1060  ELSE
1061  xi=max(xi,abs(val(j))) ! Xi = abs(max) off-diagonal element
1062  END IF
1063  END DO
1064 
1065  delta=eps*max(1.0_mpd,dgamma+xi)
1066  sn=1.0_mpd
1067  IF(n > 1) sn=1.0_mpd/sqrt(real(n*n-1,mpd))
1068  beta=sqrt(max(eps,dgamma,xi*sn)) ! beta
1069  WRITE(*,*) ' DELTA and BETA ',delta,beta
1070 
1071  DO k=2,n
1072  mk=k-ilptr(k)+ilptr(k-1)+1
1073 
1074  theta=0.0_mpd
1075 
1076  DO j=mk,k
1077  mj=j-ilptr(j)+ilptr(j-1)+1
1078  kj=ilptr(k)-k+j ! index kj
1079 
1080  DO i=max(mj,mk),j-1
1081  val(kj)=val(kj) & ! L_kj := L_kj - L_ki D_ii L_ji
1082  -val(ilptr(k)-k+i)*val(ilptr(i))*val(ilptr(j)-j+i)
1083 
1084  END DO !
1085 
1086  theta=max(theta,abs(val(kj))) ! maximum value of row
1087 
1088  IF(j /= k) THEN
1089  IF(val(ilptr(j)) /= 0.0_mpd) THEN
1090  val(kj)=val(kj)/val(ilptr(j))
1091  ELSE
1092  val(kj)=0.0_mpd
1093  END IF
1094  END IF ! L_kj := L_kj/D_jj ! D_kk
1095 
1096  IF(j == k) THEN
1097  valkj=val(kj)
1098  IF(k <= in) THEN
1099  val(kj)=max(abs(val(kj)),(theta/beta)**2,delta)
1100  IF(valkj /= val(kj)) THEN
1101  WRITE(*,*) ' Index K=',k
1102  WRITE(*,*) ' ',valkj,val(kj), (theta/beta)**2,delta,theta
1103  END IF
1104  END IF
1105  END IF
1106  END DO ! J
1107 
1108  END DO ! K
1109 
1110  DO k=1,n
1111  IF(val(ilptr(k)) /= 0.0_mpd) val(ilptr(k))=1.0_mpd/val(ilptr(k))
1112  END DO
1113 
1114  RETURN
1115 END SUBROUTINE vabdec
1116 
1122 
1123 SUBROUTINE vabmmm(n,val,ilptr)
1124  USE mpdef
1125 
1126  IMPLICIT NONE
1127  INTEGER(mpi) :: k
1128  INTEGER(mpi) :: kp
1129  INTEGER(mpi) :: kr
1130  INTEGER(mpi) :: ks
1131  INTEGER(mpi), INTENT(IN) :: n
1132  REAL(mpd), INTENT(IN OUT) :: val(*)
1133  INTEGER(mpi), INTENT(IN) :: ilptr(n)
1134  kr=1
1135  ks=1
1136  kp=1
1137  DO k=1,n
1138  IF(val(ilptr(k)) > val(ilptr(ks))) ks=k
1139  IF(val(ilptr(k)) < val(ilptr(kr))) kr=k
1140  IF(val(ilptr(k)) > 0.0.AND.val(ilptr(k)) < val(ilptr(kp))) kp=k
1141  END DO
1142  WRITE(*,*) ' Index value ',ks,val(ilptr(ks))
1143  WRITE(*,*) ' Index value ',kp,val(ilptr(kp))
1144  WRITE(*,*) ' Index value ',kr,val(ilptr(kr))
1145 
1146  RETURN
1147 END SUBROUTINE vabmmm
1148 
1159 
1160 SUBROUTINE vabslv(n,val,ilptr,x)
1161  USE mpdef
1162 
1163  IMPLICIT NONE
1164  INTEGER(mpi) :: j
1165  INTEGER(mpi) :: k
1166  INTEGER(mpi) :: mk
1167 
1168  INTEGER(mpi), INTENT(IN) :: n
1169  REAL(mpd), INTENT(IN OUT) :: val(*)
1170  INTEGER(mpi), INTENT(IN) :: ilptr(n)
1171  REAL(mpd), INTENT(IN OUT) :: x(n)
1172  ! ...
1173  DO k=2,n ! forward loop
1174  mk=k-ilptr(k)+ilptr(k-1)+1
1175  DO j=mk,k-1
1176  x(k)=x(k)-val(ilptr(k)-k+j)*x(j) ! X_k := X_k - L_kj B_j
1177  END DO
1178  END DO ! K
1179 
1180  DO k=1,n ! divide by diagonal elements
1181  x(k)=x(k)*val(ilptr(k)) ! X_k := X_k*D_kk
1182  END DO
1183 
1184  DO k=n,2,-1 ! backward loop
1185  mk=k-ilptr(k)+ilptr(k-1)+1
1186  DO j=mk,k-1
1187  x(j)=x(j)-val(ilptr(k)-k+j)*x(k) ! X_j := X_j - L_kj X_k
1188  END DO
1189  END DO ! K
1190 END SUBROUTINE vabslv
1191 
1192 ! matrix/vector products -------------------------------------------
1193 
1200 
1201 REAL(mpd) FUNCTION dbdot(n,dx,dy)
1202  USE mpdef
1203 
1204  IMPLICIT NONE
1205  INTEGER(mpi) :: i
1206  REAL(mpd) :: dtemp
1207 
1208  INTEGER(mpi), INTENT(IN) :: n
1209  REAL(mpd), INTENT(IN) :: dx(*)
1210  REAL(mpd), INTENT(IN) :: dy(*)
1211  ! ...
1212  dtemp=0.0_mpd
1213  DO i = 1,mod(n,5)
1214  dtemp=dtemp+dx(i)*dy(i)
1215  END DO
1216  DO i =mod(n,5)+1,n,5
1217  dtemp=dtemp+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2) &
1218  +dx(i+3)*dy(i+3)+dx(i+4)*dy(i+4)
1219  END DO
1220  dbdot=dtemp
1221 END FUNCTION dbdot
1222 
1231 
1232 SUBROUTINE dbaxpy(n,da,dx,dy)
1233  USE mpdef
1234 
1235  IMPLICIT NONE
1236  INTEGER(mpi) :: i
1237 
1238  INTEGER(mpi), INTENT(IN) :: n
1239  REAL(mpd), INTENT(IN) :: dx(*)
1240  REAL(mpd), INTENT(IN OUT) :: dy(*)
1241  REAL(mpd), INTENT(IN) :: da
1242  ! ...
1243  DO i=1,mod(n,4)
1244  dy(i)=dy(i)+da*dx(i)
1245  END DO
1246  DO i=mod(n,4)+1,n,4
1247  dy(i )=dy(i )+da*dx(i )
1248  dy(i+1)=dy(i+1)+da*dx(i+1)
1249  dy(i+2)=dy(i+2)+da*dx(i+2)
1250  dy(i+3)=dy(i+3)+da*dx(i+3)
1251  END DO
1252 END SUBROUTINE dbaxpy
1253 
1262 
1263 SUBROUTINE dbsvx(v,a,b,n) !
1264  USE mpdef
1265 
1266  IMPLICIT NONE
1267  INTEGER(mpi) :: i
1268  INTEGER(mpi) :: ij
1269  INTEGER(mpi) :: ijs
1270  INTEGER(mpi) :: j
1271 
1272  ! B := V * A
1273  ! N N*N N
1274 
1275  INTEGER(mpi), INTENT(IN) :: n
1276  REAL(mpd), INTENT(IN) :: v(*)
1277  REAL(mpd), INTENT(IN) :: a(*)
1278  REAL(mpd), INTENT(OUT) :: b(*)
1279 
1280  REAL(mpd) :: dsum
1281  ijs=1
1282  DO i=1,n
1283  dsum=0.0
1284  ij=ijs
1285  DO j=1,n
1286  dsum=dsum+v(ij)*a(j)
1287  IF(j < i) THEN
1288  ij=ij+1
1289  ELSE
1290  ij=ij+j
1291  END IF
1292  END DO
1293  b(i)=dsum
1294  ijs=ijs+i
1295  END DO
1296 END SUBROUTINE dbsvx
1297 
1306 
1307 SUBROUTINE dbsvxl(v,a,b,n) ! LARGE symm. matrix, vector
1308  USE mpdef
1309 
1310  IMPLICIT NONE
1311  INTEGER(mpi) :: i
1312  INTEGER(mpi) :: j
1313 
1314  ! B := V * A
1315  ! N N*N N
1316 
1317  INTEGER(mpi), INTENT(IN) :: n
1318  REAL(mpd), INTENT(IN) :: v(*)
1319  REAL(mpd), INTENT(IN) :: a(*)
1320  REAL(mpd), INTENT(OUT) :: b(*)
1321 
1322  REAL(mpd) :: dsum
1323  INTEGER(mpl) :: ij
1324  INTEGER(mpl) :: ijs
1325  ijs=1
1326  DO i=1,n
1327  dsum=0.0
1328  ij=ijs
1329  DO j=1,n
1330  dsum=dsum+v(ij)*a(j)
1331  IF(j < i) THEN
1332  ij=ij+1
1333  ELSE
1334  ij=ij+int(j,mpl)
1335  END IF
1336  END DO
1337  b(i)=dsum
1338  ijs=ijs+int(i,mpl)
1339  END DO
1340 END SUBROUTINE dbsvxl
1341 
1349 
1350 SUBROUTINE dbgax(a,x,y,m,n)
1351  USE mpdef
1352 
1353  IMPLICIT NONE
1354  INTEGER(mpi) :: i
1355  INTEGER(mpi) :: ij
1356  INTEGER(mpi) :: j
1357 
1358  REAL(mpd), INTENT(IN) :: a(*)
1359  REAL(mpd), INTENT(IN) :: x(*)
1360  REAL(mpd), INTENT(OUT) :: y(*)
1361  INTEGER(mpi), INTENT(IN) :: m
1362  INTEGER(mpi), INTENT(IN) :: n
1363 
1364  ! ...
1365  ij=0
1366  DO i=1,m
1367  y(i)=0.0_mpd
1368  DO j=1,n
1369  ij=ij+1
1370  y(i)=y(i)+a(ij)*x(j)
1371  END DO
1372  END DO
1373 END SUBROUTINE dbgax
1374 
1387 SUBROUTINE dbavat(v,a,w,n,ms)
1388  USE mpdef
1389 
1390  IMPLICIT NONE
1391  INTEGER(mpi) :: i
1392  INTEGER(mpi) :: ij
1393  INTEGER(mpi) :: ijs
1394  INTEGER(mpi) :: il
1395  INTEGER(mpi) :: j
1396  INTEGER(mpi) :: jk
1397  INTEGER(mpi) :: k
1398  INTEGER(mpi) :: l
1399  INTEGER(mpi) :: lk
1400  INTEGER(mpi) :: lkl
1401  INTEGER(mpi) :: m
1402 
1403  REAL(mpd), INTENT(IN) :: v(*)
1404  REAL(mpd), INTENT(IN) :: a(*)
1405  REAL(mpd), INTENT(INOUT) :: w(*)
1406  INTEGER(mpi), INTENT(IN) :: n
1407  INTEGER(mpi), INTENT(IN) :: ms
1408 
1409  REAL(mpd) :: cik
1410  ! ...
1411  m=ms
1412  IF (m > 0) THEN
1413  DO i=1,(m*m+m)/2
1414  w(i)=0.0_mpd ! reset output matrix
1415  END DO
1416  ELSE
1417  m=-m
1418  END IF
1419 
1420  il=-n
1421  ijs=0
1422  DO i=1,m ! do I
1423  ijs=ijs+i-1 !
1424  il=il+n !
1425  lkl=0 !
1426  DO k=1,n ! do K
1427  cik=0.0_mpd !
1428  lkl=lkl+k-1 !
1429  lk=lkl !
1430  DO l=1,k ! do L
1431  lk=lk+1 ! .
1432  cik=cik+a(il+l)*v(lk) ! .
1433  END DO ! end do L
1434  DO l=k+1,n ! do L
1435  lk=lk+l-1 ! .
1436  cik=cik+a(il+l)*v(lk) ! .
1437  END DO ! end do L
1438  jk=k !
1439  ij=ijs !
1440  DO j=1,i ! do J
1441  ij=ij+1 ! .
1442  w(ij)=w(ij)+cik*a(jk) ! .
1443  jk=jk+n ! .
1444  END DO ! end do J
1445  END DO ! end do K
1446  END DO ! end do I
1447 END SUBROUTINE dbavat
1448 
1469 SUBROUTINE dbavats(v,a,is,w,n,ms,sc)
1470  USE mpdef
1471 
1472  IMPLICIT NONE
1473  INTEGER(mpi) :: i
1474  INTEGER(mpi) :: ic
1475  INTEGER(mpi) :: ij
1476  INTEGER(mpi) :: ijs
1477  INTEGER(mpi) :: in
1478  INTEGER(mpi) :: ir
1479  INTEGER(mpi) :: j
1480  INTEGER(mpi) :: k
1481  INTEGER(mpi) :: l
1482  INTEGER(mpi) :: lk
1483  INTEGER(mpi) :: m
1484 
1485  REAL(mpd), INTENT(IN) :: v(*)
1486  REAL(mpd), INTENT(IN) :: a(*)
1487  INTEGER(mpi), INTENT(IN) :: is(*)
1488  REAL(mpd), INTENT(INOUT) :: w(*)
1489  INTEGER(mpi), INTENT(IN) :: n
1490  INTEGER(mpi), INTENT(IN) :: ms
1491  INTEGER(mpi), INTENT(OUT) :: sc(*)
1492 
1493  REAL(mpd) :: cik
1494  ! ...
1495  m=ms
1496  IF (m > 0) THEN
1497  DO i=1,(m*m+m)/2
1498  w(i)=0.0_mpd ! reset output matrix
1499  END DO
1500  ELSE
1501  m=-m
1502  END IF
1503 
1504  ! offsets in V
1505  sc(1)=0
1506  DO k=2,n
1507  sc(k)=sc(k-1)+k-1
1508  END DO
1509 
1510  ijs=0
1511  DO i=1,m
1512  ijs=ijs+i-1
1513  DO k=1,n
1514  cik=0.0_mpd
1515  DO l=is(i)+1,is(i+1),2
1516  ic=is(l)
1517  in=is(l+1)
1518  lk=sc(max(k,ic))+min(k,ic)
1519  cik=cik+a(in)*v(lk)
1520  END DO
1521  DO j=is(m+k)+1,is(m+k+1),2
1522  ir=is(j)
1523  in=is(j+1)
1524  IF (ir > i) EXIT
1525  ij=ijs+ir
1526  w(ij)=w(ij)+cik*a(in)
1527  END DO
1528  END DO
1529  END DO
1530 END SUBROUTINE dbavats
1531 
1541 
1542 SUBROUTINE dbmprv(lun,x,v,n)
1543  USE mpdef
1544 
1545  IMPLICIT NONE
1546  INTEGER(mpi) :: i
1547  INTEGER(mpi) :: ii
1548  INTEGER(mpi) :: ij
1549  INTEGER(mpi) :: j
1550  INTEGER(mpi) :: jj
1551  INTEGER(mpi) :: l
1552  INTEGER(mpi) :: m
1553  INTEGER(mpi) :: mc(15)
1554  REAL(mps) :: pd
1555  REAL(mps) :: rho
1556  REAL(mps) :: err
1557 
1558  INTEGER(mpi), INTENT(IN) :: lun
1559  REAL(mpd), INTENT(IN) :: x(*)
1560  REAL(mpd), INTENT(IN) :: v(*)
1561  INTEGER(mpi), INTENT(IN) :: n
1562 
1563  WRITE(lun,103)
1564  WRITE(lun,101)
1565  ii=0
1566  DO i=1,n
1567  ij=ii
1568  ii=ii+i
1569  err=0.0
1570  IF(v(ii) > 0.0) err=sqrt(real(v(ii),mps))
1571  l=0
1572  jj=0
1573  DO j=1,i
1574  jj=jj+j
1575  ij=ij+1
1576  rho=0.0
1577  pd=real(v(ii)*v(jj),mps)
1578  IF(pd > 0.0) rho=real(v(ij),mps)/sqrt(pd)
1579  l=l+1
1580  mc(l)=nint(100.0*abs(rho),mpi)
1581  IF(rho < 0.0) mc(l)=-mc(l)
1582  IF(j == i.OR.l == 15) THEN
1583  IF(j <= 15) THEN
1584  IF(j == i) THEN
1585  WRITE(lun,102) i,x(i),err,(mc(m),m=1,l-1)
1586  ELSE
1587  WRITE(lun,102) i,x(i),err,(mc(m),m=1,l)
1588  END IF
1589  ELSE
1590  IF(j == i) THEN
1591  WRITE(lun,103) (mc(m),m=1,l-1)
1592  ELSE
1593  WRITE(lun,103) (mc(m),m=1,l)
1594  END IF
1595  l=0
1596  END IF
1597  END IF
1598  END DO
1599  END DO
1600  WRITE(lun,104)
1601  ! 100 RETURN
1602  RETURN
1603 101 FORMAT(9x,'Param',7x,'error',7x,'correlation coefficients'/)
1604 102 FORMAT(1x,i5,2g12.4,1x,15i5)
1605 103 FORMAT(31x,15i5)
1606 104 FORMAT(33x,'(correlation coefficients in percent)')
1607 END SUBROUTINE dbmprv
1608 
1616 
1617 SUBROUTINE dbprv(lun,v,n)
1618  USE mpdef
1619 
1620  IMPLICIT NONE
1621  INTEGER(mpi), PARAMETER :: istp=6
1622  INTEGER(mpi) :: i
1623  INTEGER(mpi) :: ip
1624  INTEGER(mpi) :: ipe
1625  INTEGER(mpi) :: ipn
1626  INTEGER(mpi) :: ips
1627  INTEGER(mpi) :: k
1628 
1629  INTEGER(mpi), INTENT(IN) :: lun
1630  REAL(mpd), INTENT(IN) :: v(*)
1631  INTEGER(mpi), INTENT(IN) :: n
1632 
1633  WRITE(lun,101)
1634 
1635  DO i=1,n
1636  ips=(i*i-i)/2
1637  ipe=ips+i
1638  ip =ips
1639 100 CONTINUE
1640  ipn=ip+istp
1641  WRITE(lun,102) i, ip+1-ips, (v(k),k=ip+1,min(ipn,ipe))
1642  IF (ipn < ipe) THEN
1643  ip=ipn
1644  GO TO 100
1645  END IF
1646 END DO
1647 RETURN
1648 101 FORMAT(1x,'--- DBPRV -----------------------------------')
1649 102 FORMAT(1x,2i3,6g12.4)
1650 END SUBROUTINE dbprv
1651 
1652 ! sort -------------------------------------------------------------
1653 
1660 
1661 SUBROUTINE heapf(a,n)
1662  USE mpdef
1663 
1664  IMPLICIT NONE
1665  !
1666  INTEGER(mpi) :: i
1667  INTEGER(mpi) :: j
1668  INTEGER(mpi) :: l
1669  INTEGER(mpi) :: r
1670  REAL(mps) :: at ! pivot key value
1671 
1672  REAL(mps), INTENT(IN OUT) :: a(*)
1673  INTEGER(mpi), INTENT(IN) :: n
1674  ! ...
1675  IF(n <= 1) RETURN
1676  l=n/2+1
1677  r=n
1678 10 IF(l > 1) THEN
1679  l=l-1
1680  at =a(l)
1681  ELSE
1682  at =a(r)
1683  a(r)=a(1)
1684  r=r-1
1685  IF(r == 1) THEN
1686  a(1)=at
1687  RETURN
1688  END IF
1689  END IF
1690  i=l
1691  j=l+l
1692 20 IF(j <= r) THEN
1693  IF(j < r) THEN
1694  IF(a(j) < a(j+1)) j=j+1
1695  END IF
1696  IF(at < a(j)) THEN
1697  a(i)=a(j)
1698  i=j
1699  j=j+j
1700  ELSE
1701  j=r+1
1702  END IF
1703  GO TO 20
1704  END IF
1705  a(i)=at
1706  GO TO 10
1707 END SUBROUTINE heapf
1708 
1715 
1716 SUBROUTINE sort1k(a,n)
1717  USE mpdef
1718 
1719  IMPLICIT NONE
1720  INTEGER(mpi) :: nlev ! stack size
1721  parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1722  INTEGER(mpi) :: i
1723  INTEGER(mpi) :: j
1724  INTEGER(mpi) :: l
1725  INTEGER(mpi) :: r
1726  INTEGER(mpi) :: lev
1727  INTEGER(mpi) :: lr(nlev)
1728  INTEGER(mpi) :: lrh
1729  INTEGER(mpi) :: maxlev
1730  INTEGER(mpi) :: a1 ! pivot key
1731  INTEGER(mpi) :: at ! pivot key
1732 
1733  INTEGER(mpi), INTENT(IN OUT) :: a(*)
1734  INTEGER(mpi), INTENT(IN) :: n
1735  ! ...
1736  IF (n <= 0) RETURN
1737  maxlev=0
1738  lev=0
1739  l=1
1740  r=n
1741 10 IF(r-l == 1) THEN ! sort two elements L and R
1742  IF(a(l) > a(r)) THEN
1743  at=a(l) ! exchange L <-> R
1744  a(l)=a(r)
1745  a(r)=at
1746  END IF
1747  r=l
1748  END IF
1749  IF(r == l) THEN
1750  IF(lev <= 0) THEN
1751  ! WRITE(*,*) 'SORT1K (quicksort): maxlevel used/available =',
1752  ! + MAXLEV,'/64'
1753  RETURN
1754  END IF
1755  lev=lev-2
1756  l=lr(lev+1)
1757  r=lr(lev+2)
1758  ELSE
1759  ! LRH=(L+R)/2
1760  lrh=(l/2)+(r/2) ! avoid bit overflow
1761  IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1762  a1=a(lrh) ! middle
1763  i=l-1 ! find limits [J,I] with [L,R]
1764  j=r+1
1765 20 i=i+1
1766  IF(a(i) < a1) GO TO 20
1767 30 j=j-1
1768  IF(a(j) > a1) GO TO 30
1769  IF(i <= j) THEN
1770  at=a(i) ! exchange I <-> J
1771  a(i)=a(j)
1772  a(j)=at
1773  GO TO 20
1774  END IF
1775  IF(lev+2 > nlev) THEN
1776  CALL peend(33,'Aborted, stack overflow in quicksort')
1777  stop 'SORT1K (quicksort): stack overflow'
1778  END IF
1779  IF(r-i < j-l) THEN
1780  lr(lev+1)=l
1781  lr(lev+2)=j
1782  l=i
1783  ELSE
1784  lr(lev+1)=i
1785  lr(lev+2)=r
1786  r=j
1787  END IF
1788  lev=lev+2
1789  maxlev=max(maxlev,lev)
1790  END IF
1791  GO TO 10
1792 END SUBROUTINE sort1k
1793 
1800 
1801 SUBROUTINE sort2k(a,n)
1802  USE mpdef
1803 
1804  IMPLICIT NONE
1805  INTEGER(mpi) :: nlev ! stack size
1806  parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1807  INTEGER(mpi) :: i
1808  INTEGER(mpi) ::j
1809  INTEGER(mpi) ::l
1810  INTEGER(mpi) ::r
1811  INTEGER(mpi) ::lev
1812  INTEGER(mpi) ::lr(nlev)
1813  INTEGER(mpi) ::lrh
1814  INTEGER(mpi) ::maxlev
1815  INTEGER(mpi) ::a1 ! pivot key
1816  INTEGER(mpi) ::a2 ! pivot key
1817  INTEGER(mpi) ::at ! pivot key
1818 
1819  INTEGER(mpi), INTENT(IN OUT) :: a(2,*)
1820  INTEGER(mpi), INTENT(IN) :: n
1821  ! ...
1822  maxlev=0
1823  lev=0
1824  l=1
1825  r=n
1826 10 IF(r-l == 1) THEN ! sort two elements L and R
1827  IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
1828  at=a(1,l) ! exchange L <-> R
1829  a(1,l)=a(1,r)
1830  a(1,r)=at
1831  at=a(2,l)
1832  a(2,l)=a(2,r)
1833  a(2,r)=at
1834  END IF
1835  r=l
1836  END IF
1837  IF(r == l) THEN
1838  IF(lev <= 0) THEN
1839  WRITE(*,*) 'SORT2K (quicksort): maxlevel used/available =', maxlev,'/64'
1840  RETURN
1841  END IF
1842  lev=lev-2
1843  l=lr(lev+1)
1844  r=lr(lev+2)
1845  ELSE
1846  ! LRH=(L+R)/2
1847  lrh=(l/2)+(r/2) ! avoid bit overflow
1848  IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1849  a1=a(1,lrh) ! middle
1850  a2=a(2,lrh)
1851  i=l-1 ! find limits [J,I] with [L,R]
1852  j=r+1
1853 20 i=i+1
1854  IF(a(1,i) < a1) GO TO 20
1855  IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
1856 30 j=j-1
1857  IF(a(1,j) > a1) GO TO 30
1858  IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
1859  IF(i <= j) THEN
1860  at=a(1,i) ! exchange I <-> J
1861  a(1,i)=a(1,j)
1862  a(1,j)=at
1863  at=a(2,i)
1864  a(2,i)=a(2,j)
1865  a(2,j)=at
1866  GO TO 20
1867  END IF
1868  IF(lev+2 > nlev) THEN
1869  CALL peend(33,'Aborted, stack overflow in quicksort')
1870  stop 'SORT2K (quicksort): stack overflow'
1871  END IF
1872  IF(r-i < j-l) THEN
1873  lr(lev+1)=l
1874  lr(lev+2)=j
1875  l=i
1876  ELSE
1877  lr(lev+1)=i
1878  lr(lev+2)=r
1879  r=j
1880  END IF
1881  lev=lev+2
1882  maxlev=max(maxlev,lev)
1883  END IF
1884  GO TO 10
1885 END SUBROUTINE sort2k
1886 
1893 
1894 SUBROUTINE sort2i(a,n)
1895  USE mpdef
1896 
1897  IMPLICIT NONE
1898  INTEGER(mpi) :: nlev ! stack size
1899  parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1900  INTEGER(mpi) :: i
1901  INTEGER(mpi) ::j
1902  INTEGER(mpi) ::l
1903  INTEGER(mpi) ::r
1904  INTEGER(mpi) ::lev
1905  INTEGER(mpi) ::lr(nlev)
1906  INTEGER(mpi) ::lrh
1907  INTEGER(mpi) ::maxlev
1908  INTEGER(mpi) ::a1 ! pivot key
1909  INTEGER(mpi) ::a2 ! pivot key
1910  INTEGER(mpi) ::at(3)
1911 
1912  INTEGER(mpi), INTENT(IN OUT) :: a(3,*)
1913  INTEGER(mpi), INTENT(IN) :: n
1914  ! ...
1915  maxlev=0
1916  lev=0
1917  l=1
1918  r=n
1919 10 IF(r-l == 1) THEN ! sort two elements L and R
1920  IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
1921  at=a(:,l) ! exchange L <-> R
1922  a(:,l)=a(:,r)
1923  a(:,r)=at
1924  END IF
1925  r=l
1926  END IF
1927  IF(r == l) THEN
1928  IF(lev <= 0) THEN
1929  WRITE(*,*) 'SORT2I (quicksort): maxlevel used/available =', maxlev,'/64'
1930  RETURN
1931  END IF
1932  lev=lev-2
1933  l=lr(lev+1)
1934  r=lr(lev+2)
1935  ELSE
1936  ! LRH=(L+R)/2
1937  lrh=(l/2)+(r/2) ! avoid bit overflow
1938  IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1939  a1=a(1,lrh) ! middle
1940  a2=a(2,lrh)
1941  i=l-1 ! find limits [J,I] with [L,R]
1942  j=r+1
1943 20 i=i+1
1944  IF(a(1,i) < a1) GO TO 20
1945  IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
1946 30 j=j-1
1947  IF(a(1,j) > a1) GO TO 30
1948  IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
1949  IF(i <= j) THEN
1950  IF(a(1,i) == a(1,j).AND.a(2,i) == a(2,j)) GO TO 20 ! equal -> keep order
1951  at=a(:,i) ! exchange I <-> J
1952  a(:,i)=a(:,j)
1953  a(:,j)=at
1954  GO TO 20
1955  END IF
1956  IF(lev+2 > nlev) THEN
1957  CALL peend(33,'Aborted, stack overflow in quicksort')
1958  stop 'SORT2I (quicksort): stack overflow'
1959  END IF
1960  IF(r-i < j-l) THEN
1961  lr(lev+1)=l
1962  lr(lev+2)=j
1963  l=i
1964  ELSE
1965  lr(lev+1)=i
1966  lr(lev+2)=r
1967  r=j
1968  END IF
1969  lev=lev+2
1970  maxlev=max(maxlev,lev)
1971  END IF
1972  GO TO 10
1973 END SUBROUTINE sort2i
1974 
1981 
1982 SUBROUTINE sort22(a,n)
1983  USE mpdef
1984 
1985  IMPLICIT NONE
1986  INTEGER(mpi) :: nlev ! stack size
1987  parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1988  INTEGER(mpi) :: i
1989  INTEGER(mpi) ::j
1990  INTEGER(mpi) ::l
1991  INTEGER(mpi) ::r
1992  INTEGER(mpi) ::lev
1993  INTEGER(mpi) ::lr(nlev)
1994  INTEGER(mpi) ::lrh
1995  INTEGER(mpi) ::maxlev
1996  INTEGER(mpi) ::a1 ! pivot key
1997  INTEGER(mpi) ::a2 ! pivot key
1998  INTEGER(mpi) ::at(4)
1999 
2000  INTEGER(mpi), INTENT(IN OUT) :: a(4,*)
2001  INTEGER(mpi), INTENT(IN) :: n
2002  ! ...
2003  maxlev=0
2004  lev=0
2005  l=1
2006  r=n
2007  IF (n<=1) RETURN
2008 10 IF(r-l == 1) THEN ! sort two elements L and R
2009  IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
2010  at=a(:,l) ! exchange L <-> R
2011  a(:,l)=a(:,r)
2012  a(:,r)=at
2013  END IF
2014  r=l
2015  END IF
2016  IF(r == l) THEN
2017  IF(lev <= 0) THEN
2018  WRITE(*,*) 'SORT22 (quicksort): maxlevel used/available =', maxlev,'/64'
2019  RETURN
2020  END IF
2021  lev=lev-2
2022  l=lr(lev+1)
2023  r=lr(lev+2)
2024  ELSE
2025  ! LRH=(L+R)/2
2026  lrh=(l/2)+(r/2) ! avoid bit overflow
2027  IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
2028  a1=a(1,lrh) ! middle
2029  a2=a(2,lrh)
2030  i=l-1 ! find limits [J,I] with [L,R]
2031  j=r+1
2032 20 i=i+1
2033  IF(a(1,i) < a1) GO TO 20
2034  IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
2035 30 j=j-1
2036  IF(a(1,j) > a1) GO TO 30
2037  IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
2038  IF(i <= j) THEN
2039  at=a(:,i) ! exchange I <-> J
2040  a(:,i)=a(:,j)
2041  a(:,j)=at
2042  GO TO 20
2043  END IF
2044  IF(lev+2 > nlev) THEN
2045  CALL peend(33,'Aborted, stack overflow in quicksort')
2046  stop 'SORT22 (quicksort): stack overflow'
2047  END IF
2048  IF(r-i < j-l) THEN
2049  lr(lev+1)=l
2050  lr(lev+2)=j
2051  l=i
2052  ELSE
2053  lr(lev+1)=i
2054  lr(lev+2)=r
2055  r=j
2056  END IF
2057  lev=lev+2
2058  maxlev=max(maxlev,lev)
2059  END IF
2060  GO TO 10
2061 END SUBROUTINE sort22
2062 
2070 
2071 REAL(mps) FUNCTION chindl(n,nd)
2072  USE mpdef
2073 
2074  IMPLICIT NONE
2075  INTEGER(mpi) :: m
2076  INTEGER(mpi), INTENT(IN) :: n
2077  INTEGER(mpi), INTENT(IN) :: nd
2078  !
2079  REAL(mps) :: sn(3)
2080  REAL(mps) ::table(30,3)
2081  ! REAL PN(3)
2082  ! DATA PN/0.31731,0.0455002785,2.69985E-3/ ! probabilities
2083  DATA sn/0.47523,1.690140,2.782170/
2084  DATA table/ 1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, &
2085  1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, &
2086  1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, &
2087  1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040, &
2088  4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, &
2089  1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, &
2090  1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, &
2091  1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742, &
2092  9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, &
2093  2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, &
2094  2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, &
2095  1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681/
2096  SAVE sn,table
2097  ! ...
2098  IF(nd < 1) THEN
2099  chindl=0.0
2100  ELSE
2101  m=max(1,min(n,3)) ! 1, 2 or 3 sigmas
2102  IF(nd <= 30) THEN
2103  chindl=table(nd,m) ! from table
2104  ELSE ! approximation for ND > 30
2105  chindl=(sn(m)+sqrt(real(nd+nd-1,mps)))**2/real(nd+nd,mps)
2106  END IF
2107  END IF
2108 END FUNCTION chindl
2109 
2147 
2148 SUBROUTINE lltdec(n,c,india,nrkd,iopt)
2149  USE mpdef
2150 
2151  IMPLICIT NONE
2152  INTEGER(mpi) :: i
2153  INTEGER(mpi) :: j
2154  INTEGER(mpi) :: k
2155  INTEGER(mpi) :: kj
2156  INTEGER(mpi) :: mj
2157  INTEGER(mpi) :: mk
2158  REAL(mpd) ::diag
2159 
2160  INTEGER(mpi), INTENT(IN) :: n
2161  REAL(mpd), INTENT(IN OUT) :: c(*)
2162  INTEGER(mpi), INTENT(IN) :: india(n)
2163  INTEGER(mpi), INTENT(OUT) :: nrkd
2164  INTEGER(mpi), INTENT(IN) :: iopt
2165  REAL(mpd) eps
2166  ! ...
2167  eps = 16.0_mpd * epsilon(eps) ! 16 * precision(mpd)
2168 
2169  ! ..
2170  nrkd=0
2171  diag=0.0_mpd
2172  IF(c(india(1)) > 0.0) THEN
2173  c(india(1))=1.0_mpd/sqrt(c(india(1))) ! square root
2174  ELSE
2175  c(india(1))=0.0_mpd
2176  nrkd=-1
2177  END IF
2178 
2179  DO k=2,n
2180  mk=k-india(k)+india(k-1)+1 ! first index in row K
2181  DO j=mk,k ! loop over row K with index J
2182  mj=1
2183  IF(j > 1) mj=j-india(j)+india(j-1)+1 ! first index in row J
2184  kj=india(k)-k+j ! index kj
2185  diag=c(india(j)) ! j-th diagonal element
2186 
2187  DO i=max(mj,mk),j-1
2188  ! L_kj = L_kj - L_ki *D_ii *L_ji
2189  c(kj)=c(kj) - c(india(k)-k+i)*c(india(j)-j+i)
2190  END DO ! I
2191 
2192  IF(j /= k) c(kj)=c(kj)*diag
2193  END DO ! J
2194 
2195  IF(c(india(k)) > eps*diag) THEN ! test for linear dependence
2196  c(india(k))=1.0_mpd/sqrt(c(india(k))) ! square root
2197  ELSE
2198  DO j=mk,k ! reset row K
2199  c(india(k)-k+j)=0.0_mpd
2200  END DO ! J
2201  IF (iopt > 0 .and. diag > 0.0) THEN ! skyline
2202  c(india(k))=1.0_mpd/sqrt(diag) ! square root
2203  ELSE
2204  IF(nrkd == 0) THEN
2205  nrkd=-k
2206  ELSE
2207  IF(nrkd < 0) nrkd=1
2208  nrkd=nrkd+1
2209  END IF
2210  END IF
2211  END IF
2212 
2213  END DO ! K
2214  RETURN
2215 END SUBROUTINE lltdec
2216 
2217 
2229 
2230 SUBROUTINE lltfwd(n,c,india,x)
2231  USE mpdef
2232 
2233  IMPLICIT NONE
2234  INTEGER(mpi) :: j
2235  INTEGER(mpi) :: k
2236 
2237  INTEGER(mpi), INTENT(IN) :: n
2238  REAL(mpd), INTENT(IN) :: c(*)
2239  INTEGER(mpi), INTENT(IN) :: india(n)
2240  REAL(mpd), INTENT(IN OUT) :: x(n)
2241 
2242  x(1)=x(1)*c(india(1))
2243  DO k=2,n ! forward loop
2244  DO j=k-india(k)+india(k-1)+1,k-1
2245  x(k)=x(k)-c(india(k)-k+j)*x(j) ! X_k := X_k - L_kj * B_j
2246  END DO ! J
2247  x(k)=x(k)*c(india(k))
2248  END DO ! K
2249 
2250  RETURN
2251 END SUBROUTINE lltfwd
2252 
2264 
2265 SUBROUTINE lltbwd(n,c,india,x)
2266  USE mpdef
2267 
2268  IMPLICIT NONE
2269  INTEGER(mpi) :: j
2270  INTEGER(mpi) :: k
2271 
2272  INTEGER(mpi), INTENT(IN) :: n
2273  REAL(mpd), INTENT(IN) :: c(*)
2274  INTEGER(mpi), INTENT(IN) :: india(n)
2275  REAL(mpd), INTENT(IN OUT) :: x(n)
2276 
2277  DO k=n,2,-1 ! backward loop
2278  x(k)=x(k)*c(india(k))
2279  DO j=k-india(k)+india(k-1)+1,k-1
2280  x(j)=x(j)-c(india(k)-k+j)*x(k) ! X_j := X_j - L_kj * X_k
2281  END DO ! J
2282  END DO ! K
2283  x(1)=x(1)*c(india(1))
2284 
2285  RETURN
2286 END SUBROUTINE lltbwd
2287 
2305 SUBROUTINE equdec(n,m,ls,c,india,nrkd,nrkd2)
2306  USE mpdef
2307 
2308  IMPLICIT NONE
2309  INTEGER(mpi) :: i
2310  INTEGER(mpi) :: j
2311  INTEGER(mpi) :: jk
2312  INTEGER(mpi) :: k
2313 
2314  INTEGER(mpi), INTENT(IN) :: n
2315  INTEGER(mpi), INTENT(IN) :: m
2316  INTEGER(mpi), INTENT(IN) :: ls
2317  REAL(mpd), INTENT(IN OUT) :: c(*)
2318  INTEGER(mpi), INTENT(IN OUT) :: india(n+m)
2319  INTEGER(mpi), INTENT(OUT) :: nrkd
2320  INTEGER(mpi), INTENT(OUT) :: nrkd2
2321 
2322  ! ...
2323 
2324  nrkd=0
2325  nrkd2=0
2326 
2327  CALL lltdec(n,c,india,nrkd,ls) ! decomposition G G^T
2328 
2329  IF (m>0) THEN
2330  DO i=1,m
2331  CALL lltfwd(n,c,india,c(india(n)+(i-1)*n+1)) ! forward solution K
2332  END DO
2333 
2334  jk=india(n)+n*m
2335  DO j=1,m
2336  DO k=1,j
2337  jk=jk+1
2338  c(jk)=0.0_mpd ! product K K^T
2339  DO i=1,n
2340  c(jk)=c(jk)+c(india(n)+(j-1)*n+i)*c(india(n)+(k-1)*n+i)
2341  END DO
2342  END DO
2343  END DO
2344 
2345  india(n+1)=1
2346  DO i=2,m
2347  india(n+i)=india(n+i-1)+min(i,m) ! pointer for K K^T
2348  END DO
2349 
2350  CALL lltdec(m,c(india(n)+n*m+1),india(n+1),nrkd2,0) ! decomp. H H^T
2351  ENDIF
2352 
2353  RETURN
2354 END SUBROUTINE equdec
2355 
2371 SUBROUTINE equslv(n,m,c,india,x) ! solution vector
2372  USE mpdef
2373 
2374  IMPLICIT NONE
2375  INTEGER(mpi) :: i
2376  INTEGER(mpi) :: j
2377 
2378  INTEGER(mpi), INTENT(IN) :: n
2379  INTEGER(mpi), INTENT(IN) :: m
2380  REAL(mpd), INTENT(IN) :: c(*)
2381  INTEGER(mpi), INTENT(IN) :: india(n+m)
2382  REAL(mpd), INTENT(IN OUT) :: x(n+m)
2383 
2384  CALL lltfwd(n,c,india,x) ! result is u
2385 
2386  IF (m>0) THEN
2387  DO i=1,m
2388  DO j=1,n
2389  x(n+i)=x(n+i)-x(j)*c(india(n)+(i-1)*n+j) ! g - K u
2390  END DO
2391  END DO
2392  CALL lltfwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is v
2393 
2394 
2395  CALL lltbwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is -y
2396  DO i=1,m
2397  x(n+i)=-x(n+i) ! result is +y
2398  END DO
2399 
2400  DO i=1,n
2401  DO j=1,m
2402  x(i)=x(i)-x(n+j)*c(india(n)+(j-1)*n+i) ! u - K^T y
2403  END DO
2404  END DO
2405  ENDIF
2406 
2407  CALL lltbwd(n,c,india,x) ! result is x
2408 
2409  RETURN
2410 END SUBROUTINE equslv
2411 
2441 
2442 SUBROUTINE precon(p,n,c,cu,a,s,nrkd)
2443  USE mpdef
2444 
2445  IMPLICIT NONE
2446  INTEGER(mpi) :: i
2447  INTEGER(mpi) :: ii
2448  INTEGER(mpi) :: j
2449  INTEGER(mpi) :: jj
2450  INTEGER(mpi) :: jk
2451  INTEGER(mpi) :: k
2452  INTEGER(mpi) :: kk
2453 
2454  INTEGER(mpi), INTENT(IN) :: p
2455  INTEGER(mpi), INTENT(IN) :: n
2456  REAL(mpd), INTENT(IN) :: c(n)
2457  REAL(mpd), INTENT(OUT) :: cu(n)
2458  REAL(mpd), INTENT(IN OUT) :: a(n,p)
2459  REAL(mpd), INTENT(OUT) :: s((p*p+p)/2)
2460  INTEGER(mpi), INTENT(OUT) :: nrkd
2461 
2462  REAL(mpd) :: div
2463  REAL(mpd) :: ratio
2464 
2465  nrkd=0
2466  DO i=1,(p*p+p)/2
2467  s(i)=0.0_mpd
2468  END DO
2469  DO i=1,n
2470  jk=0
2471  div=c(i) ! copy
2472  IF (div > 0.0_mpd) THEN
2473  cu(i)=1.0_mpd/sqrt(div)
2474  ELSE
2475  cu(i)=0.0_mpd
2476  nrkd=nrkd+1
2477  END IF
2478  DO j=1,p
2479  a(i,j)=a(i,j)*cu(i) ! K = A C^{-1/2}
2480  DO k=1,j
2481  jk=jk+1
2482  s(jk)=s(jk)+a(i,j)*a(i,k) ! S = symmetric matrix K K^T
2483  END DO
2484  END DO
2485  END DO
2486 
2487  ii=0
2488  DO i=1,p ! S -> H D H^T (Cholesky)
2489  ii=ii+i
2490  IF(s(ii) /= 0.0_mpd) s(ii)=1.0_mpd/s(ii)
2491  jj=ii
2492  DO j=i+1,p
2493  ratio=s(i+jj)*s(ii)
2494  kk=jj
2495  DO k=j,p
2496  s(kk+j)=s(kk+j)-s(kk+i)*ratio
2497  kk=kk+k
2498  END DO ! K
2499  s(i+jj)=ratio
2500  jj=jj+j
2501  END DO ! J
2502  END DO ! I
2503  RETURN
2504 END SUBROUTINE precon
2505 
2515 
2516 SUBROUTINE presol(p,n,cu,a,s,x,y) ! solution
2517  USE mpdef
2518 
2519  IMPLICIT NONE
2520  INTEGER(mpi) :: i
2521  INTEGER(mpi) :: j
2522  INTEGER(mpi) :: jj
2523  INTEGER(mpi) :: k
2524  INTEGER(mpi) :: kk
2525 
2526  INTEGER(mpi), INTENT(IN) :: p
2527  INTEGER(mpi), INTENT(IN) :: n
2528 
2529  REAL(mpd), INTENT(IN) :: cu(n)
2530  REAL(mpd), INTENT(IN) :: a(n,p)
2531  REAL(mpd), INTENT(IN) :: s((p*p+p)/2)
2532  REAL(mpd), INTENT(OUT) :: x(n+p)
2533  REAL(mpd), INTENT(IN) :: y(n+p)
2534 
2535  REAL(mpd) :: dsum
2536 
2537  DO i=1,n+p
2538  x(i)=y(i)
2539  END DO
2540  DO i=1,n
2541  x(i)=x(i)*cu(i) ! u =C^{-1/2} y
2542  DO j=1,p
2543  x(n+j)=x(n+j)-a(i,j)*x(i) ! d - K u
2544  END DO
2545  END DO
2546 
2547  jj=0
2548  DO j=1,p ! Cholesky solution for v
2549  dsum=x(n+j)
2550  DO k=1,j-1
2551  dsum=dsum-s(k+jj)*x(n+k) ! H v = d - K u
2552  END DO
2553  x(n+j)=dsum ! -> v
2554  jj=jj+j
2555  END DO
2556 
2557  DO j=p,1,-1 ! solution for lambda
2558  dsum=x(n+j)*s(jj)
2559  kk=jj
2560  DO k=j+1,p
2561  dsum=dsum+s(kk+j)*x(n+k) ! D H^T lambda = -v
2562  kk=kk+k
2563  END DO
2564  x(n+j)=-dsum ! -> lambda
2565  jj=jj-j
2566  END DO
2567 
2568  DO i=1,n ! u - K^T lambda
2569  DO j=1,p
2570  x(i)=x(i)-a(i,j)*x(n+j)
2571  END DO
2572  END DO
2573  DO i=1,n
2574  x(i)=x(i)*cu(i) ! x = C^{-1/2} u
2575  END DO
2576 
2577 END SUBROUTINE presol
2578 
2579 
2580 ! 090817 C. Kleinwort, DESY-FH1
2625 SUBROUTINE sqmibb(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
2626  USE mpdef
2627 
2628  ! REAL(mpd) scratch arrays:
2629  ! VBND(N*(NBND+1)) = storage of band part
2630  ! VBDR(N* NBDR) = storage of border part
2631  ! AUX (N* NBDR) = intermediate results
2632 
2633  ! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only)
2634 
2635  IMPLICIT NONE
2636  INTEGER(mpi) :: i
2637  INTEGER(mpi) :: ib
2638  INTEGER(mpi) :: ij
2639  INTEGER(mpi) :: ioff
2640  INTEGER(mpi) :: ip
2641  INTEGER(mpi) :: ip1
2642  INTEGER(mpi) :: ip2
2643  INTEGER(mpi) :: is
2644  INTEGER(mpi) :: j
2645  INTEGER(mpi) :: j0
2646  INTEGER(mpi) :: jb
2647  INTEGER(mpi) :: joff
2648  INTEGER(mpi) :: mp1
2649  INTEGER(mpi) :: nb1
2650  INTEGER(mpi) :: nmb
2651  INTEGER(mpi) :: npri
2652  INTEGER(mpi) :: nrankb
2653 
2654  REAL(mpd), INTENT(IN OUT) :: v(*)
2655  REAL(mpd), INTENT(OUT) :: b(n)
2656  INTEGER(mpi), INTENT(IN) :: n
2657  INTEGER(mpi), INTENT(IN) :: nbdr
2658  INTEGER(mpi), INTENT(IN) :: nbnd
2659  INTEGER(mpi), INTENT(IN) :: inv
2660  INTEGER(mpi), INTENT(OUT) :: nrank
2661 
2662  REAL(mpd), INTENT(OUT) :: vbnd(n*(nbnd+1))
2663  REAL(mpd), INTENT(OUT) :: vbdr(n*nbdr)
2664  REAL(mpd), INTENT(OUT) :: aux(n*nbdr)
2665  REAL(mpd), INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
2666  REAL(mpd), INTENT(OUT) :: vzru(nbdr)
2667  REAL(mpd), INTENT(OUT) :: scdiag(nbdr)
2668  INTEGER(mpi), INTENT(OUT) :: scflag(nbdr)
2669 
2670  SAVE npri
2671  DATA npri / 100 /
2672  ! ...
2673  nrank=0
2674  nb1=nbdr+1
2675  mp1=nbnd+1
2676  nmb=n-nbdr
2677  ! copy band part
2678  DO i=nb1,n
2679  ip=(i*(i+1))/2
2680  is=0
2681  DO j=i,min(n,i+nbnd)
2682  ip=ip+is
2683  is=j
2684  ib=j-i+1
2685  vbnd(ib+(i-nb1)*mp1)=v(ip)
2686  END DO
2687  END DO
2688  ! copy border part
2689  IF (nbdr > 0) THEN
2690  ioff=0
2691  DO i=1,nbdr
2692  ip=(i*(i+1))/2
2693  is=0
2694  DO j=i,n
2695  ip=ip+is
2696  is=j
2697  vbdr(ioff+j)=v(ip)
2698  END DO
2699  ioff=ioff+n
2700  END DO
2701  END IF
2702 
2703  CALL dbcdec(vbnd,mp1,nmb,aux)
2704  ! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable
2705  ! CALL DBCPRB(VBND,MP1,NMB)
2706  ip=1
2707  DO i=1, nmb
2708  IF (vbnd(ip) <= 0.0_mpd) THEN
2709  npri=npri-1
2710  IF (npri >= 0) THEN
2711  IF (vbnd(ip) == 0.0_mpd) THEN
2712  print *, ' SQMIBB matrix singular', n, nbdr, nbnd
2713  ELSE
2714  print *, ' SQMIBB matrix not positive definite', n, nbdr, nbnd
2715  END IF
2716  END IF
2717  ! return zeros
2718  DO ip=1,n
2719  b(ip)=0.0_mpd
2720  END DO
2721  DO ip=1,(n*n+n)/2
2722  v(ip)=0.0_mpd
2723  END DO
2724  RETURN
2725  END IF
2726  ip=ip+mp1
2727  END DO
2728  nrank=nmb
2729 
2730  IF (nbdr == 0) THEN ! special case NBDR=0
2731 
2732  CALL dbcslv(vbnd,mp1,nmb,b,b)
2733  IF (inv > 0) THEN
2734  IF (inv > 1) THEN
2735  CALL dbcinv(vbnd,mp1,nmb,v)
2736  ELSE
2737  CALL dbcinb(vbnd,mp1,nmb,v)
2738  END IF
2739  END IF
2740 
2741  ELSE ! general case NBDR>0
2742 
2743  ioff=nb1
2744  DO ib=1,nbdr
2745  ! solve for aux. vectors
2746  CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff),aux(ioff))
2747  ! zT ru
2748  vzru(ib)=b(ib)
2749  DO i=0,nmb-1
2750  vzru(ib)=vzru(ib)-b(nb1+i)*aux(ioff+i)
2751  END DO
2752  ioff=ioff+n
2753  END DO
2754  ! solve for band part only
2755  CALL dbcslv(vbnd,mp1,nmb,b(nb1),b(nb1))
2756  ! Ck - cT z
2757  ip=0
2758  ioff=nb1
2759  DO ib=1,nbdr
2760  joff=nb1
2761  DO jb=1,ib
2762  ip=ip+1
2763  vbk(ip)=v(ip)
2764  DO i=0,nmb-1
2765  vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
2766  END DO
2767  joff=joff+n
2768  END DO
2769  ioff=ioff+n
2770  END DO
2771  ! solve border part
2772  CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
2773  IF (nrankb == nbdr) THEN
2774  nrank=nrank+nbdr
2775  ELSE
2776  npri=npri-1
2777  IF (npri >= 0) print *, ' SQMIBB undef border ', n, nbdr, nbnd, nrankb
2778  DO ib=1,nbdr
2779  vzru(ib)=0.0_mpd
2780  END DO
2781  DO ip=(nbdr*nbdr+nbdr)/2,1,-1
2782  vbk(ip)=0.0_mpd
2783  END DO
2784  END IF
2785  ! smoothed data points
2786  ioff=nb1
2787  DO ib=1, nbdr
2788  b(ib) = vzru(ib)
2789  DO i=0,nmb-1
2790  b(nb1+i)=b(nb1+i)-b(ib)*aux(ioff+i)
2791  END DO
2792  ioff=ioff+n
2793  END DO
2794  ! inverse requested ?
2795  IF (inv > 0) THEN
2796  IF (inv > 1) THEN
2797  CALL dbcinv(vbnd,mp1,nmb,v)
2798  ELSE
2799  CALL dbcinb(vbnd,mp1,nmb,v)
2800  END IF
2801  ! expand/correct from NMB to N
2802  ip1=(nmb*nmb+nmb)/2
2803  ip2=(n*n+n)/2
2804  DO i=nmb-1,0,-1
2805  j0=0
2806  IF (inv == 1) j0=max(0,i-nbnd)
2807  DO j=i,j0,-1
2808  v(ip2)=v(ip1)
2809  ioff=nb1
2810  DO ib=1,nbdr
2811  joff=nb1
2812  DO jb=1,nbdr
2813  ij=max(ib,jb)
2814  ij=(ij*ij-ij)/2+min(ib,jb)
2815  v(ip2)=v(ip2)+vbk(ij)*aux(ioff+i)*aux(joff+j)
2816  joff=joff+n
2817  END DO
2818  ioff=ioff+n
2819  END DO
2820  ip1=ip1-1
2821  ip2=ip2-1
2822  END DO
2823  ip1=ip1-j0
2824  ip2=ip2-j0
2825 
2826  DO ib=nbdr,1,-1
2827  v(ip2)=0.0_mpd
2828  joff=nb1
2829  DO jb=1,nbdr
2830  ij=max(ib,jb)
2831  ij=(ij*ij-ij)/2+min(ib,jb)
2832  v(ip2)=v(ip2)-vbk(ij)*aux(i+joff)
2833  joff=joff+n
2834  END DO
2835  ip2=ip2-1
2836  END DO
2837  END DO
2838 
2839  DO ip=(nbdr*nbdr+nbdr)/2,1,-1
2840  v(ip2)=vbk(ip)
2841  ip2=ip2-1
2842  END DO
2843 
2844  END IF
2845  END IF
2846 
2847 END SUBROUTINE sqmibb
2848 
2849 ! 181105 C. Kleinwort, DESY-BELLE
2881 SUBROUTINE sqmibb2(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
2882  USE mpdef
2883 
2884  ! REAL(mpd) scratch arrays:
2885  ! VBND(N*(NBND+1)) = storage of band part
2886  ! VBDR(N* NBDR) = storage of border part
2887  ! AUX (N* NBDR) = intermediate results
2888 
2889  ! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only)
2890 
2891  IMPLICIT NONE
2892  INTEGER(mpi) :: i
2893  INTEGER(mpi) :: ib
2894  INTEGER(mpi) :: ij
2895  INTEGER(mpi) :: ioff
2896  INTEGER(mpi) :: ip
2897  INTEGER(mpi) :: ip1
2898  INTEGER(mpi) :: is
2899  INTEGER(mpi) :: j
2900  INTEGER(mpi) :: j0
2901  INTEGER(mpi) :: jb
2902  INTEGER(mpi) :: joff
2903  INTEGER(mpi) :: koff
2904  INTEGER(mpi) :: mp1
2905  INTEGER(mpi) :: nb1
2906  INTEGER(mpi) :: nmb
2907  INTEGER(mpi) :: npri
2908  INTEGER(mpi) :: nrankb
2909 
2910  REAL(mpd), INTENT(IN OUT) :: v(*)
2911  REAL(mpd), INTENT(OUT) :: b(n)
2912  INTEGER(mpi), INTENT(IN) :: n
2913  INTEGER(mpi), INTENT(IN) :: nbdr
2914  INTEGER(mpi), INTENT(IN) :: nbnd
2915  INTEGER(mpi), INTENT(IN) :: inv
2916  INTEGER(mpi), INTENT(OUT) :: nrank
2917 
2918  REAL(mpd), INTENT(OUT) :: vbnd(n*(nbnd+1))
2919  REAL(mpd), INTENT(OUT) :: vbdr(n*nbdr)
2920  REAL(mpd), INTENT(OUT) :: aux(n*nbdr)
2921  REAL(mpd), INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
2922  REAL(mpd), INTENT(OUT) :: vzru(nbdr)
2923  REAL(mpd), INTENT(OUT) :: scdiag(nbdr)
2924  INTEGER(mpi), INTENT(OUT) :: scflag(nbdr)
2925 
2926  SAVE npri
2927  DATA npri / 100 /
2928  ! ...
2929  nrank=0
2930  mp1=nbnd+1
2931  nmb=n-nbdr
2932  nb1=nmb+1
2933  ! copy band part
2934  DO i=1,nmb
2935  ip=(i*(i+1))/2
2936  is=0
2937  DO j=i,min(nmb,i+nbnd)
2938  ip=ip+is
2939  is=j
2940  ib=j-i+1
2941  vbnd(ib+(i-1)*mp1)=v(ip)
2942  END DO
2943  END DO
2944  ! copy border part
2945  IF (nbdr > 0) THEN
2946  ioff=0
2947  DO i=nb1,n
2948  ip=(i*(i-1))/2
2949  DO j=1,i
2950  vbdr(ioff+j)=v(ip+j)
2951  END DO
2952  ioff=ioff+n
2953  END DO
2954  END IF
2955 
2956  CALL dbcdec(vbnd,mp1,nmb,aux)
2957  ! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable
2958  ! CALL DBCPRB(VBND,MP1,NMB)
2959  ip=1
2960  DO i=1, nmb
2961  IF (vbnd(ip) <= 0.0_mpd) THEN
2962  npri=npri-1
2963  IF (npri >= 0) THEN
2964  IF (vbnd(ip) == 0.0_mpd) THEN
2965  print *, ' SQMIBB2 matrix singular', n, nbdr, nbnd
2966  ELSE
2967  print *, ' SQMIBB2 matrix not positive definite', n, nbdr, nbnd
2968  END IF
2969  END IF
2970  ! return zeros
2971  DO ip=1,n
2972  b(ip)=0.0_mpd
2973  END DO
2974  DO ip=1,(n*n+n)/2
2975  v(ip)=0.0_mpd
2976  END DO
2977  RETURN
2978  END IF
2979  ip=ip+mp1
2980  END DO
2981  nrank=nmb
2982 
2983  IF (nbdr == 0) THEN ! special case NBDR=0
2984 
2985  CALL dbcslv(vbnd,mp1,nmb,b,b)
2986  IF (inv > 0) THEN
2987  IF (inv > 1) THEN
2988  CALL dbcinv(vbnd,mp1,nmb,v)
2989  ELSE
2990  CALL dbcinb(vbnd,mp1,nmb,v)
2991  END IF
2992  END IF
2993 
2994  ELSE ! general case NBDR>0
2995 
2996  ioff=0
2997  DO ib=1,nbdr
2998  ! solve for aux. vectors
2999  CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff+1),aux(ioff+1))
3000  ! zT ru
3001  vzru(ib)=b(nmb+ib)
3002  DO i=1,nmb
3003  vzru(ib)=vzru(ib)-b(i)*aux(ioff+i)
3004  END DO
3005  ioff=ioff+n
3006  END DO
3007  ! solve for band part only
3008  CALL dbcslv(vbnd,mp1,nmb,b,b)
3009  ! Ck - cT z
3010  ip=0
3011  ioff=0
3012  koff=nmb
3013  DO ib=1,nbdr
3014  joff=0
3015  DO jb=1,ib
3016  ip=ip+1
3017  vbk(ip)=vbdr(koff+jb)
3018  DO i=1,nmb
3019  vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
3020  END DO
3021  joff=joff+n
3022  END DO
3023  ioff=ioff+n
3024  koff=koff+n
3025  END DO
3026 
3027  ! solve border part
3028  CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
3029  IF (nrankb == nbdr) THEN
3030  nrank=nrank+nbdr
3031  ELSE
3032  npri=npri-1
3033  IF (npri >= 0) print *, ' SQMIBB2 undef border ', n, nbdr, nbnd, nrankb
3034  DO ib=1,nbdr
3035  vzru(ib)=0.0_mpd
3036  END DO
3037  DO ip=(nbdr*nbdr+nbdr)/2,1,-1
3038  vbk(ip)=0.0_mpd
3039  END DO
3040  END IF
3041  ! smoothed data points
3042  ioff=0
3043  DO ib=1, nbdr
3044  DO i=1,nmb
3045  b(i)=b(i)-vzru(ib)*aux(ioff+i)
3046  END DO
3047  ioff=ioff+n
3048  b(nmb+ib)=vzru(ib)
3049  END DO
3050  ! inverse requested ?
3051  IF (inv > 0) THEN
3052  IF (inv > 1) THEN
3053  CALL dbcinv(vbnd,mp1,nmb,v)
3054  ELSE
3055  CALL dbcinb(vbnd,mp1,nmb,v)
3056  END IF
3057  ! assemble band and border
3058  IF (nbdr > 0) THEN
3059  ! band part
3060  ip1=(nmb*nmb+nmb)/2
3061  DO i=nmb-1,0,-1
3062  j0=0
3063  IF (inv == 1) j0=max(0,i-nbnd)
3064  DO j=i,j0,-1
3065  ioff=1
3066  DO ib=1,nbdr
3067  joff=1
3068  DO jb=1,nbdr
3069  ij=max(ib,jb)
3070  ij=(ij*ij-ij)/2+min(ib,jb)
3071  v(ip1)=v(ip1)+vbk(ij)*aux(ioff+i)*aux(joff+j)
3072  joff=joff+n
3073  END DO
3074  ioff=ioff+n
3075  END DO
3076  ip1=ip1-1
3077  END DO
3078  ip1=ip1-j0
3079  END DO
3080  ! border part
3081  ip1=(nmb*nmb+nmb)/2
3082  ip=0
3083  DO ib=1,nbdr
3084  DO i=1,nmb
3085  ip1=ip1+1
3086  v(ip1)=0.0_mpd
3087  joff=0
3088  DO jb=1,nbdr
3089  ij=max(ib,jb)
3090  ij=(ij*ij-ij)/2+min(ib,jb)
3091  v(ip1)=v(ip1)-vbk(ij)*aux(i+joff)
3092  joff=joff+n
3093  END DO
3094  END DO
3095  DO jb=1,ib
3096  ip1=ip1+1
3097  ip=ip+1
3098  v(ip1)=vbk(ip)
3099  END DO
3100  END DO
3101 
3102  END IF
3103  END IF
3104  END IF
3105 
3106 END SUBROUTINE sqmibb2
subroutine dbcdec(w, mp1, n, aux)
Decomposition of symmetric band matrix.
subroutine monpgs(i)
Progress monitoring.
Definition: mpmon.f90:68
subroutine cholin(g, v, n)
Inversion after decomposition.
Definition: mpnum.f90:837
subroutine sqmibb2(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Band bordered matrix.
Definition: mpnum.f90:2882
subroutine vabmmm(n, val, ilptr)
Variable band matrix print minimum and maximum.
Definition: mpnum.f90:1124
subroutine dbavats(v, a, is, w, n, ms, sc)
A V AT product (similarity, sparse).
Definition: mpnum.f90:1470
subroutine sqminl(v, b, n, nrank, diag, next, vk, mon)
Matrix inversion for LARGE matrices.
Definition: mpnum.f90:230
subroutine precon(p, n, c, cu, a, s, nrkd)
Constrained preconditioner, decomposition.
Definition: mpnum.f90:2443
subroutine devsol(n, diag, u, b, x, work)
Solution by diagonalization.
Definition: mpnum.f90:649
subroutine dbsvxl(v, a, b, n)
Product LARGE symmetric matrix, vector.
Definition: mpnum.f90:1308
subroutine devrot(n, diag, u, v, work, iwork)
Diagonalization.
Definition: mpnum.f90:369
subroutine dbmprv(lun, x, v, n)
Print symmetric matrix, vector.
Definition: mpnum.f90:1543
subroutine dbaxpy(n, da, dx, dy)
Multiply, addition.
Definition: mpnum.f90:1233
subroutine equslv(n, m, c, india, x)
Solution of equilibrium systems (after decomposition).
Definition: mpnum.f90:2372
subroutine heapf(a, n)
Heap sort direct (real).
Definition: mpnum.f90:1662
real(mpd) function dbdot(n, dx, dy)
Dot product.
Definition: mpnum.f90:1202
subroutine equdec(n, m, ls, c, india, nrkd, nrkd2)
Decomposition of equilibrium systems.
Definition: mpnum.f90:2306
subroutine vabdec(n, val, ilptr)
Variable band matrix decomposition.
Definition: mpnum.f90:1012
subroutine chslv2(g, x, n)
Solve A*x=b using Cholesky decomposition.
Definition: mpnum.f90:953
subroutine sort22(a, n)
Quick sort 2 with index.
Definition: mpnum.f90:1983
subroutine vabslv(n, val, ilptr, x)
Variable band matrix solution.
Definition: mpnum.f90:1161
subroutine choldc(g, n)
Cholesky decomposition.
Definition: mpnum.f90:745
subroutine sort1k(a, n)
Quick sort 1.
Definition: mpnum.f90:1717
subroutine sqminv(v, b, n, nrank, diag, next)
Matrix inversion and solution.
Definition: mpnum.f90:97
subroutine lltbwd(n, c, india, x)
Backward solution.
Definition: mpnum.f90:2266
subroutine dbgax(a, x, y, m, n)
Multiply general M-by-N matrix A and N-vector X.
Definition: mpnum.f90:1351
subroutine sqmibb(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Bordered band matrix.
Definition: mpnum.f90:2626
subroutine presol(p, n, cu, a, s, x, y)
Constrained preconditioner, solution.
Definition: mpnum.f90:2517
subroutine dbprv(lun, v, n)
Print symmetric matrix.
Definition: mpnum.f90:1618
subroutine devinv(n, diag, u, v)
Inversion by diagonalization.
Definition: mpnum.f90:696
subroutine lltdec(n, c, india, nrkd, iopt)
LLT decomposition.
Definition: mpnum.f90:2149
subroutine cholsl(g, x, n)
Solution after decomposition.
Definition: mpnum.f90:791
subroutine chdec2(g, n, nrank, evmax, evmin, mon)
Cholesky decomposition (LARGE pos.
Definition: mpnum.f90:891
subroutine sort2k(a, n)
Quick sort 2.
Definition: mpnum.f90:1802
subroutine devsig(n, diag, u, b, coef)
Calculate significances.
Definition: mpnum.f90:611
real(mps) function chindl(n, nd)
Chi2/ndf cuts.
Definition: mpnum.f90:2072
subroutine lltfwd(n, c, india, x)
Forward solution.
Definition: mpnum.f90:2231
subroutine dbsvx(v, a, b, n)
Product symmetric matrix, vector.
Definition: mpnum.f90:1264
subroutine dbavat(v, a, w, n, ms)
A V AT product (similarity).
Definition: mpnum.f90:1388
subroutine sort2i(a, n)
Quick sort 2 with index.
Definition: mpnum.f90:1895
Definition of constants.
Definition: mpdef.f90:24
integer, parameter mpd
Definition: mpdef.f90:38
integer, parameter mpi
Definition: mpdef.f90:34
integer, parameter mpl
Definition: mpdef.f90:36
integer, parameter parameter
Definition: mpdef.f90:34
integer, parameter mps
Definition: mpdef.f90:37
subroutine peend(icode, cmessage)
Print exit code.
Definition: pede.f90:10469