Millepede-II V04-16-01
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
76
77!----------------------------------------------------------------------
96
97SUBROUTINE sqminv(v,b,n,nrank,diag,next) ! matrix inversion
98 USE mpdef
99
100 IMPLICIT NONE
101 INTEGER(mpi) :: i
102 INTEGER(mpi) :: ij
103 INTEGER(mpi) :: j
104 INTEGER(mpi) :: jj
105 INTEGER(mpi) :: jk
106 INTEGER(mpi) :: jl
107 INTEGER(mpi) :: k
108 INTEGER(mpi) :: kk
109 INTEGER(mpi) :: l
110 INTEGER(mpi) :: last
111 INTEGER(mpi) :: lk
112 INTEGER(mpi) :: next0
113
114 INTEGER(mpi), INTENT(IN) :: n
115 REAL(mpd), INTENT(IN OUT) :: v((n*n+n)/2)
116 REAL(mpd), INTENT(OUT) :: b(n)
117 INTEGER(mpi), INTENT(OUT) :: nrank
118 REAL(mpd), INTENT(OUT) :: diag(n)
119 INTEGER(mpi), INTENT(OUT) :: next(n)
120 REAL(mpd) :: vkk
121 REAL(mpd) :: vjk
122
123 !REAL(mpd), PARAMETER :: eps=1.0E-10_mpd
124 REAL(mpd) eps
125 ! ...
126 eps = 16.0_mpd * epsilon(eps) ! 16 * precision(mpd)
127
128 next0=1
129 l=1
130 DO i=1,n
131 next(i)=i+1 ! set "next" pointer
132 diag(i)=abs(v((i*i+i)/2)) ! save abs of diagonal elements
133 END DO
134 next(n)=-1 ! end flag
135
136 nrank=0
137 DO i=1,n ! start of loop
138 k =0
139 vkk=0.0_mpd
140
141 j=next0
142 last=0
14305 IF(j > 0) THEN
144 jj=(j*j+j)/2
145 IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
146 vkk=v(jj)
147 k=j
148 l=last
149 END IF
150 last=j
151 j=next(last)
152 GO TO 05
153 END IF
154
155 IF(k /= 0) THEN ! pivot found
156 kk=(k*k+k)/2
157 IF(l == 0) THEN
158 next0=next(k)
159 ELSE
160 next(l)=next(k)
161 END IF
162 next(k)=0 ! index is used, reset
163 nrank=nrank+1 ! increase rank and ...
164 vkk =1.0/vkk
165 v(kk) =-vkk
166 b(k) =b(k)*vkk
167 jk =kk-k
168 jl =0
169 DO j=1,n ! elimination
170 IF(j == k) THEN
171 jk=kk
172 jl=jl+j
173 ELSE
174 IF(j < k) THEN
175 jk=jk+1
176 ELSE
177 jk=jk+j-1
178 END IF
179 vjk =v(jk)
180 v(jk)=vkk*vjk
181 b(j) =b(j)-b(k)*vjk
182 lk =kk-k
183 DO l=1,j
184 jl=jl+1
185 IF(l == k) THEN
186 lk=kk
187 ELSE
188 IF(l < k) THEN
189 lk=lk+1
190 ELSE
191 lk=lk+l-1
192 END IF
193 v(jl)=v(jl)-v(lk)*vjk
194 END IF
195 END DO
196 END IF
197 END DO
198 ELSE
199 DO k=1,n
200 IF(next(k) /= 0) THEN
201 b(k)=0.0_mpd ! clear vector element
202 DO j=1,k
203 IF(next(j) /= 0) v((k*k-k)/2+j)=0.0_mpd ! clear matrix row/col
204 END DO
205 END IF
206 END DO
207 GO TO 10
208 END IF
209 END DO ! end of loop
210 10 DO ij=1,(n*n+n)/2
211 v(ij)=-v(ij) ! finally reverse sign of all matrix elements
212 END DO
213END SUBROUTINE sqminv
214
229
230SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk,mon) !
231 USE mpdef
232
233 IMPLICIT NONE
234 INTEGER(mpi) :: i
235 INTEGER(mpi) :: j
236 INTEGER(mpi) :: k
237 INTEGER(mpi) :: l
238 INTEGER(mpi) :: last
239 INTEGER(mpi) :: next0
240
241 INTEGER(mpi), INTENT(IN) :: n
242 REAL(mpd), INTENT(IN OUT) :: v((INT(n,mpl)*INT(n,mpl)+INT(n,mpl))/2)
243 REAL(mpd), INTENT(OUT) :: b(n)
244 INTEGER(mpi), INTENT(OUT) :: nrank
245 REAL(mpd), INTENT(OUT) :: diag(n)
246 INTEGER(mpi), INTENT(OUT) :: next(n)
247 REAL(mpd), INTENT(OUT) :: vk(n)
248 INTEGER(mpi), INTENT(IN) :: mon
249
250 INTEGER(mpl) :: i8
251 INTEGER(mpl) :: j8
252 INTEGER(mpl) :: jj
253 INTEGER(mpl) :: k8
254 INTEGER(mpl) :: kk
255 INTEGER(mpl) :: kkmk
256 INTEGER(mpl) :: jk
257 INTEGER(mpl) :: jl
258
259 REAL(mpd) :: vkk
260 REAL(mpd) :: vjk
261
262 REAL(mpd), PARAMETER :: eps=1.0e-10_mpd
263 ! ...
264 next0=1
265 l=1
266 DO i=1,n
267 i8=int(i,mpl)
268 next(i)=i+1 ! set "next" pointer
269 diag(i)=abs(v((i8*i8+i8)/2)) ! save abs of diagonal elements
270 END DO
271 next(n)=-1 ! end flag
272
273 nrank=0
274 DO i=1,n ! start of loop
275 ! monitoring ?
276 IF(mon>0) CALL monpgs(i)
277 k =0
278 vkk=0.0_mpd
279 j=next0
280 last=0
28105 IF(j > 0) THEN
282 j8=int(j,mpl)
283 jj=(j8*j8+j8)/2
284 IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
285 vkk=v(jj)
286 k=j
287 l=last
288 END IF
289 last=j
290 j=next(last)
291 GO TO 05
292 END IF
293
294 IF(k /= 0) THEN ! pivot found
295 k8=int(k,mpl)
296 kk=(k8*k8+k8)/2
297 kkmk=kk-k8
298 IF(l == 0) THEN
299 next0=next(k)
300 ELSE
301 next(l)=next(k)
302 END IF
303 next(k)=0 ! index is used, reset
304 nrank=nrank+1 ! increase rank and ...
305 vkk =1.0/vkk
306 v(kk) =-vkk
307 b(k) =b(k)*vkk
308 ! elimination
309 jk =kkmk
310 DO j=1,n
311 IF(j == k) THEN
312 jk=kk
313 vk(j)=0.
314 ELSE
315 IF(j < k) THEN
316 jk=jk+1
317 ELSE
318 jk=jk+int(j,mpl)-1
319 END IF
320 v(jk)=v(jk)*vkk
321 vk(j)=v(jk)
322 END IF
323 END DO
324 ! parallelize row loop
325 ! slot of 128 'J' for next idle thread (optimized on Intel Xeon)
326 !$OMP PARALLEL DO &
327 !$OMP PRIVATE(JL,VJK,J8) &
328 !$OMP SCHEDULE(DYNAMIC,128)
329 DO j=n,1,-1
330 IF(j == k) cycle
331 j8=int(j,mpl)
332 jl=j8*(j8-1)/2
333 vjk =vk(j)/vkk
334 b(j) =b(j)-b(k)*vjk
335 v(jl+1:jl+j)=v(jl+1:jl+j)-vk(1:j)*vjk
336 END DO
337 !$OMP END PARALLEL DO
338 ELSE
339 DO k=1,n
340 k8=int(k,mpl)
341 kk=(k8*k8-k8)/2
342 IF(next(k) /= 0) THEN
343 b(k)=0.0_mpd ! clear vector element
344 DO j=1,k
345 IF(next(j) /= 0) v(kk+int(j,mpl))=0.0_mpd ! clear matrix row/col
346 END DO
347 END IF
348 END DO
349 GO TO 10
350 END IF
351 END DO ! end of loop
352 10 DO jj=1,(int(n,mpl)*int(n+1,mpl))/2
353 v(jj)=-v(jj) ! finally reverse sign of all matrix elements
354 END DO
355END SUBROUTINE sqminl
356
368
369SUBROUTINE devrot(n,diag,u,v,work,iwork) ! diagonalization
370 USE mpdef
371
372 IMPLICIT NONE
373
374 INTEGER(mpi), INTENT(IN) :: n
375 REAL(mpd), INTENT(OUT) :: diag(n)
376 REAL(mpd), INTENT(OUT) :: u(n,n)
377 REAL(mpd), INTENT(IN) :: v((n*n+n)/2)
378 REAL(mpd), INTENT(OUT) :: work(n)
379 INTEGER(mpi), INTENT(OUT) :: iwork(n)
380
381
382 INTEGER(mpi), PARAMETER :: itmax=30
383 REAL(mpd), PARAMETER :: tol=epsilon(tol)
384 REAL(mpd), PARAMETER :: eps=epsilon(eps)
385
386 REAL(mpd) :: f
387 REAL(mpd) :: g
388 REAL(mpd) :: h
389 REAL(mpd) :: sh
390 REAL(mpd) :: hh
391 REAL(mpd) :: b
392 REAL(mpd) :: p
393 REAL(mpd) :: r
394 REAL(mpd) :: s
395 REAL(mpd) :: c
396 REAL(mpd) :: workd
397
398 INTEGER(mpi) :: ij
399 INTEGER(mpi) :: i
400 INTEGER(mpi) :: j
401 INTEGER(mpi) :: k
402 INTEGER(mpi) :: l
403 INTEGER(mpi) :: m
404 INTEGER(mpi) :: ll
405 ! ...
406 ! 1. part: symmetric matrix V reduced to tridiagonal from
407 ij=0
408 DO i=1,n
409 DO j=1,i
410 ij=ij+1
411 u(i,j)=v(ij) ! copy half of symmetric matirx
412 END DO
413 END DO
414
415 DO i=n,2,-1
416 l=i-2
417 f=u(i,i-1)
418 g=0.0_mpd
419 IF(l /= 0) THEN
420 DO k=1,l
421 IF(abs(u(i,k)) > tol) g=g+u(i,k)*u(i,k)
422 END DO
423 h=g+f*f
424 END IF
425 IF(g < tol) THEN ! G too small
426 work(i)=f ! skip transformation
427 h =0.0_mpd
428 ELSE
429 l=l+1
430 sh=sqrt(h)
431 IF(f >= 0.0_mpd) sh=-sh
432 g=sh
433 work(i)=sh
434 h=h-f*g
435 u(i,i-1)=f-g
436 f=0.0_mpd
437 DO j=1,l
438 u(j,i)=u(i,j)/h
439 g=0.0_mpd
440 ! form element of a u
441 DO k=1,j
442 IF(abs(u(j,k)) > tol.AND.abs(u(i,k)) > tol) THEN
443 g=g+u(j,k)*u(i,k)
444 END IF
445 END DO
446 DO k=j+1,l
447 IF(abs(u(k,j)) > tol.AND.abs(u(i,k)) > tol) THEN
448 g=g+u(k,j)*u(i,k)
449 END IF
450 END DO
451 work(j)=g/h
452 f=f+g*u(j,i)
453 END DO
454 ! form k
455 hh=f/(h+h)
456 ! form reduced a
457 DO j=1,l
458 f=u(i,j)
459 work(j)=work(j)-hh*f
460 g=work(j)
461 DO k=1,j
462 u(j,k)=u(j,k)-f*work(k)-g*u(i,k)
463 END DO
464 END DO
465 END IF
466 diag(i)=h
467 END DO
468
469 diag(1)=0.0_mpd
470 work(1)=0.0_mpd
471
472 ! accumulation of transformation matrices
473 DO i=1,n
474 IF(diag(i) /= 0.0) THEN
475 DO j=1,i-1
476 g=0.0_mpd
477 DO k=1,i-1
478 g=g+u(i,k)*u(k,j)
479 END DO
480 DO k=1,i-1
481 u(k,j)=u(k,j)-g*u(k,i)
482 END DO
483 END DO
484 END IF
485 diag(i)=u(i,i)
486 u(i,i)=1.0_mpd
487 DO j=1,i-1
488 u(i,j)=0.0_mpd
489 u(j,i)=0.0_mpd
490 END DO
491 END DO
492
493 ! 2. part: diagonalization of tridiagonal matrix
494 DO i=2,n
495 work(i-1)=work(i)
496 END DO
497 work(n)=0.0_mpd
498 b=0.0_mpd
499 f=0.0_mpd
500
501 DO l=1,n
502 j=0
503 h=eps*(abs(diag(l))+abs(work(l)))
504 IF(b < h) b=h
505 DO m=l,n
506 IF(abs(work(m)) <= b) GO TO 10 ! look for small sub-diagonal element
507 END DO
508 m=l
50910 IF(m == l) GO TO 30
510 ! next iteration
51120 IF(j == itmax) THEN
512 WRITE(*,*) 'DEVROT: Iteration limit reached'
513 CALL peend(32,'Aborted, iteration limit reached in diagonalization')
514 stop
515 END IF
516 j=j+1
517 g=diag(l)
518 p=(diag(l+1)-g)/(2.0_mpd*work(l))
519 r=sqrt(1.0_mpd+p*p)
520 diag(l)=work(l)
521 IF(p < 0.0_mpd) diag(l)=diag(l)/(p-r)
522 IF(p >= 0.0_mpd) diag(l)=diag(l)/(p+r)
523 h=g-diag(l)
524 DO i=l+1,n
525 diag(i)=diag(i)-h
526 END DO
527 f=f+h
528 ! QL transformation
529 p=diag(m)
530 c=1.0_mpd
531 s=0.0_mpd
532 DO i=m-1,l,-1 ! reverse loop
533 g=c*work(i)
534 h=c*p
535 IF(abs(p) >= abs(work(i))) THEN
536 c=work(i)/p
537 r=sqrt(1.0_mpd+c*c)
538 work(i+1)=s*p*r
539 s=c/r
540 c=1.0_mpd/r
541 ELSE
542 c=p/work(i)
543 r=sqrt(1.0_mpd+c*c)
544 work(i+1)=s*work(i)*r
545 s=1.0_mpd/r
546 c=c/r
547 END IF
548 p=c*diag(i)-s*g
549 diag(i+1)=h+s*(c*g+s*diag(i))
550 ! form vector
551 DO k=1,n
552 h=u(k,i+1)
553 u(k,i+1)=s*u(k,i)+c*h
554 u(k,i)=c*u(k,i)-s*h
555 END DO
556 END DO
557 work(l)=s*p
558 diag(l)=c*p
559 IF(abs(work(l)) > b) GO TO 20 ! next iteration
56030 diag(l)=diag(l)+f
561 END DO
562 DO i=1,n
563 iwork(i)=i
564 END DO
565
566 m=1
56740 m=1+3*m ! determine initial increment
568 IF(m <= n) GO TO 40
56950 m=m/3
570 DO j=1,n-m ! sort with increment M
571 l=j
57260 IF(diag(iwork(l+m)) > diag(iwork(l))) THEN ! compare
573 ll=iwork(l+m) ! exchange the two index values
574 iwork(l+m)=iwork(l)
575 iwork(l)=ll
576 l=l-m
577 IF(l > 0) GO TO 60
578 END IF
579 END DO
580 IF(m > 1) GO TO 50
581
582 DO i=1,n
583 IF(iwork(i) /= i) THEN
584 ! move vector from position I to the work area
585 workd=diag(i)
586 DO l=1,n
587 work(l)=u(l,i)
588 END DO
589 k=i
59070 j=k
591 k=iwork(j)
592 iwork(j)=j
593 IF(k /= i) THEN
594 ! move vector from position K to the (free) position J
595 diag(j)=diag(k)
596 DO l=1,n
597 u(l,j)=u(l,k)
598 END DO
599 GO TO 70
600 END IF
601 ! move vector from the work area to position J
602 diag(j)=workd
603 DO l=1,n
604 u(l,j)=work(l)
605 END DO
606 END IF
607 END DO
608END SUBROUTINE devrot
609
611SUBROUTINE devsig(n,diag,u,b,coef)
612 USE mpdef
613
614 IMPLICIT NONE
615
616 INTEGER(mpi), INTENT(IN) :: n
617 REAL(mpd), INTENT(IN) :: diag(n)
618 REAL(mpd), INTENT(IN) :: u(n,n)
619 REAL(mpd), INTENT(IN) :: b(n)
620 REAL(mpd), INTENT(OUT) :: coef(n)
621 INTEGER(mpi) :: i
622 INTEGER(mpi) :: j
623 REAL(mpd) :: dsum
624 ! ...
625 DO i=1,n
626 coef(i)=0.0_mpd
627 IF(diag(i) > 0.0_mpd) THEN
628 dsum=0.0_mpd
629 DO j=1,n
630 dsum=dsum+u(j,i)*b(j)
631 END DO
632 coef(i)=abs(dsum)/sqrt(diag(i))
633 END IF
634 END DO
635END SUBROUTINE devsig
636
637
648
649SUBROUTINE devsol(n,diag,u,b,x,work)
650 USE mpdef
651
652 IMPLICIT NONE
653
654 INTEGER(mpi), INTENT(IN) :: n
655 REAL(mpd), INTENT(IN) :: diag(n)
656 REAL(mpd), INTENT(IN) :: u(n,n)
657 REAL(mpd), INTENT(IN) :: b(n)
658 REAL(mpd), INTENT(OUT) :: x(n)
659 REAL(mpd), INTENT(OUT) :: work(n)
660 INTEGER(mpi) :: i
661 INTEGER(mpi) :: j
662 INTEGER(mpi) :: jj
663 REAL(mpd) :: s
664 ! ...
665 DO j=1,n
666 s=0.0_mpd
667 work(j)=0.0_mpd
668 IF(diag(j) /= 0.0_mpd) THEN
669 DO i=1,n
670 ! j-th eigenvector is U(.,J)
671 s=s+u(i,j)*b(i)
672 END DO
673 work(j)=s/diag(j)
674 END IF
675 END DO
676
677 DO j=1,n
678 s=0.0_mpd
679 DO jj=1,n
680 s=s+u(j,jj)*work(jj)
681 END DO
682 x(j)=s
683 END DO
684! WRITE(*,*) 'DEVSOL'
685! WRITE(*,*) 'X ',X
686END SUBROUTINE devsol
687
695
696SUBROUTINE devinv(n,diag,u,v)
697 USE mpdef
698
699 IMPLICIT NONE
700 INTEGER(mpi) :: i
701 INTEGER(mpi) :: ij
702 INTEGER(mpi) :: j
703 INTEGER(mpi) :: k
704
705 INTEGER(mpi), INTENT(IN) :: n
706 REAL(mpd), INTENT(IN) :: diag(n)
707 REAL(mpd), INTENT(IN) :: u(n,n)
708 REAL(mpd), INTENT(OUT) :: v((n*n+n)/2)
709 REAL(mpd) :: dsum
710 ! ...
711 ij=0
712 DO i=1,n
713 DO j=1,i
714 ij=ij+1
715 dsum=0.0_mpd
716 DO k=1,n
717 IF(diag(k) /= 0.0_mpd) THEN
718 dsum=dsum+u(i,k)*u(j,k)/diag(k)
719 END IF
720 END DO
721 v(ij)=dsum
722 END DO
723 END DO
724END SUBROUTINE devinv
725
726
745SUBROUTINE choldc(g,n)
746 USE mpdef
747
748 IMPLICIT NONE
749 INTEGER(mpi) :: i
750 INTEGER(mpi) :: ii
751 INTEGER(mpi) :: j
752 INTEGER(mpi) :: jj
753 INTEGER(mpi) :: k
754 INTEGER(mpi) :: kk
755
756 INTEGER(mpi), INTENT(IN) :: n
757 REAL(mpd), INTENT(IN OUT) :: g((n*n+n)/2)
758
759 REAL(mpd) :: ratio
760 ! ...
761 ii=0
762 DO i=1,n
763 ii=ii+i
764 IF(g(ii) /= 0.0) g(ii)=1.0/g(ii) ! (I,I) div !
765 jj=ii
766 DO j=i+1,n
767 ratio=g(i+jj)*g(ii) ! (I,J) (I,I)
768 kk=jj
769 DO k=j,n
770 g(kk+j)=g(kk+j)-g(kk+i)*ratio ! (K,J) (K,I)
771 kk=kk+k
772 END DO ! K
773 g(i+jj)=ratio ! (I,J)
774 jj=jj+j
775 END DO ! J
776 END DO ! I
777 RETURN
778END SUBROUTINE choldc
779
791SUBROUTINE cholsl(g,x,n)
792 USE mpdef
793
794 IMPLICIT NONE
795 REAL(mpd) :: dsum
796 INTEGER(mpi) :: i
797 INTEGER(mpi) :: ii
798 INTEGER(mpi) :: k
799 INTEGER(mpi) :: kk
800
801 INTEGER(mpi), INTENT(IN) :: n
802 REAL(mpd), INTENT(IN) :: g((n*n+n)/2)
803 REAL(mpd), INTENT(IN OUT) :: x(n)
804
805 ii=0
806 DO i=1,n
807 dsum=x(i)
808 DO k=1,i-1
809 dsum=dsum-g(k+ii)*x(k) ! (K,I)
810 END DO
811 x(i)=dsum
812 ii=ii+i
813 END DO
814 DO i=n,1,-1
815 dsum=x(i)*g(ii) ! (I,I)
816 kk=ii
817 DO k=i+1,n
818 dsum=dsum-g(kk+i)*x(k) ! (K,I)
819 kk=kk+k
820 END DO
821 x(i)=dsum
822 ii=ii-i
823 END DO
824 RETURN
825END SUBROUTINE cholsl
826
837SUBROUTINE cholin(g,v,n)
838 USE mpdef
839
840 IMPLICIT NONE
841 REAL(mpd) :: dsum
842 INTEGER(mpi) :: i
843 INTEGER(mpi) :: ii
844 INTEGER(mpi) :: j
845 INTEGER(mpi) :: k
846 INTEGER(mpi) :: l
847 INTEGER(mpi) :: m
848
849 INTEGER(mpi), INTENT(IN) :: n
850 REAL(mpd), INTENT(IN) :: g((n*n+n)/2)
851 REAL(mpd), INTENT( OUT) :: v((n*n+n)/2)
852
853 ii=(n*n-n)/2
854 DO i=n,1,-1
855 dsum=g(ii+i) ! (I,I)
856 DO j=i,1,-1
857 DO k=j+1,n
858 l=min(i,k)
859 m=max(i,k)
860 dsum=dsum-g(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K)
861 END DO
862 v(ii+j)=dsum ! (I,J)
863 dsum=0.0_mpd
864 END DO
865 ii=ii-i+1
866 END DO
867END SUBROUTINE cholin
868
869! 201026 C. Kleinwort, DESY-BELLE
891SUBROUTINE chdec2(g,n,nrank,evmax,evmin,mon)
892 USE mpdef
893
894 IMPLICIT NONE
895 INTEGER(mpi) :: i
896 INTEGER(mpi) :: j
897
898 INTEGER(mpl) :: ii
899 INTEGER(mpl) :: jj
900 REAL(mpd) :: ratio
901
902 INTEGER(mpi), INTENT(IN) :: n
903 REAL(mpd), INTENT(IN OUT) :: g((INT(n,mpl)*INT(n,mpl)+INT(n,mpl))/2)
904 INTEGER(mpi), INTENT(OUT) :: nrank
905 REAL(mpd), INTENT(OUT) :: evmin
906 REAL(mpd), INTENT(OUT) :: evmax
907 INTEGER(mpi), INTENT(IN) :: mon
908
909 nrank=0
910 ii=(int(n,mpl)*int(n+1,mpl))/2
911 DO i=n,1,-1
912 ! monitoring ?
913 IF(mon>0) CALL monpgs(n+1-i)
914 IF (g(ii) > 0.0_mpd) THEN
915 ! update rank, min, max eigenvalue
916 nrank=nrank+1
917 IF (nrank == 1) THEN
918 evmax=g(ii)
919 evmin=g(ii)
920 ELSE
921 evmax=max(evmax,g(ii))
922 evmin=min(evmin,g(ii))
923 END IF
924 g(ii)=1.0/g(ii) ! (I,I) div !
925 END IF
926 ii=ii-i
927 ! parallelize row loop
928 ! slot of 32 'J' for next idle thread (optimized on Intel Xeon)
929 !$OMP PARALLEL DO &
930 !$OMP PRIVATE(RATIO,JJ) &
931 !$OMP SCHEDULE(DYNAMIC,32)
932 DO j=1,i-1
933 ratio=g(ii+j)*g(ii+i) ! (I,J) (I,I)
934 IF (ratio == 0.0_mpd) cycle
935 jj=(int(j-1,mpl)*int(j,mpl))/2
936 g(jj+1:jj+j)=g(jj+1:jj+j)-g(ii+1:ii+j)*ratio ! (K,J) (K,I)
937 END DO ! J
938 !$OMP END PARALLEL DO
939 g(ii+1:ii+i-1)=g(ii+1:ii+i-1)*g(ii+i) ! (I,J)
940 END DO ! I
941
942END SUBROUTINE chdec2
943
944! 201026 C. Kleinwort, DESY-BELLE
953SUBROUTINE chslv2(g,x,n)
954 USE mpdef
955
956 IMPLICIT NONE
957 INTEGER(mpi) :: i
958 INTEGER(mpi) :: k
959 INTEGER(mpl) :: ii
960 INTEGER(mpl) :: kk
961 REAL(mpd) :: dsum
962
963 INTEGER(mpi), INTENT(IN) :: n
964 REAL(mpd), INTENT(IN) :: g((INT(n,mpl)*INT(n,mpl)+INT(n,mpl))/2)
965 REAL(mpd), INTENT(IN OUT) :: x(n)
966
967 ii=(int(n,mpl)*int(n+1,mpl))/2
968 DO i=n,1,-1
969 dsum=x(i)
970 kk=ii
971 DO k=i+1,n
972 dsum=dsum-g(kk+i)*x(k) ! (K,I)
973 kk=kk+k
974 END DO
975 x(i)=dsum
976 ii=ii-i
977 END DO
978 DO i=1,n
979 dsum=x(i)*g(ii+i) ! (I,I)
980 DO k=1,i-1
981 dsum=dsum-g(k+ii)*x(k) ! (K,I)
982 END DO
983 x(i)=dsum
984 ii=ii+i
985 END DO
986
987END SUBROUTINE chslv2
988
989! variable band matrix operations ----------------------------------
990
1011
1012SUBROUTINE vabdec(n,val,ilptr)
1013 USE mpdef
1014
1015 IMPLICIT NONE
1016
1017 INTEGER(mpi) :: i
1018 INTEGER(mpi) :: in
1019 INTEGER(mpi) :: j
1020 INTEGER(mpi) :: k
1021 INTEGER(mpi) :: kj
1022 INTEGER(mpi) :: mj
1023 INTEGER(mpi) :: mk
1024 REAL(mpd) :: sn
1025 REAL(mpd) :: beta
1026 REAL(mpd) :: delta
1027 REAL(mpd) :: theta
1028
1029 INTEGER(mpi), INTENT(IN) :: n
1030 REAL(mpd), INTENT(IN OUT) :: val(*)
1031 INTEGER(mpi), INTENT(IN) :: ilptr(n)
1032
1033 REAL(mpd) :: dgamma
1034 REAL(mpd) :: xi
1035 REAL(mpd) :: valkj
1036
1037 REAL(mpd), PARAMETER :: one=1.0_mpd
1038 REAL(mpd), PARAMETER :: two=2.0_mpd
1039 REAL(mpd), PARAMETER :: eps = epsilon(eps)
1040
1041 WRITE(*,*) 'Variable band matrix Cholesky decomposition'
1042
1043 dgamma=0.0_mpd
1044 i=1
1045 DO j=1,ilptr(n) ! loop thrugh all matrix elements
1046 IF(ilptr(i) == j) THEN ! diagonal element
1047 IF(val(j) <= 0.0_mpd) GO TO 01 ! exit loop for negative diag
1048 dgamma=max(dgamma,abs(val(j))) ! max diagonal element
1049 i=i+1
1050 END IF
1051 END DO
1052 i=n+1
105301 in=i-1 ! IN positive diagonal elements
1054 WRITE(*,*) ' ',in,' positive diagonal elements'
1055 xi=0.0_mpd
1056 i=1
1057 DO j=1,ilptr(in) ! loop for positive diagonal elements
1058 ! through all matrix elements
1059 IF(ilptr(i) == j) THEN ! diagonal element
1060 i=i+1
1061 ELSE
1062 xi=max(xi,abs(val(j))) ! Xi = abs(max) off-diagonal element
1063 END IF
1064 END DO
1065
1066 delta=eps*max(1.0_mpd,dgamma+xi)
1067 sn=1.0_mpd
1068 IF(n > 1) sn=1.0_mpd/sqrt(real(n*n-1,mpd))
1069 beta=sqrt(max(eps,dgamma,xi*sn)) ! beta
1070 WRITE(*,*) ' DELTA and BETA ',delta,beta
1071
1072 DO k=2,n
1073 mk=k-ilptr(k)+ilptr(k-1)+1
1074
1075 theta=0.0_mpd
1076
1077 DO j=mk,k
1078 mj=j-ilptr(j)+ilptr(j-1)+1
1079 kj=ilptr(k)-k+j ! index kj
1080
1081 DO i=max(mj,mk),j-1
1082 val(kj)=val(kj) & ! L_kj := L_kj - L_ki D_ii L_ji
1083 -val(ilptr(k)-k+i)*val(ilptr(i))*val(ilptr(j)-j+i)
1084
1085 END DO !
1086
1087 theta=max(theta,abs(val(kj))) ! maximum value of row
1088
1089 IF(j /= k) THEN
1090 IF(val(ilptr(j)) /= 0.0_mpd) THEN
1091 val(kj)=val(kj)/val(ilptr(j))
1092 ELSE
1093 val(kj)=0.0_mpd
1094 END IF
1095 END IF ! L_kj := L_kj/D_jj ! D_kk
1096
1097 IF(j == k) THEN
1098 valkj=val(kj)
1099 IF(k <= in) THEN
1100 val(kj)=max(abs(val(kj)),(theta/beta)**2,delta)
1101 IF(valkj /= val(kj)) THEN
1102 WRITE(*,*) ' Index K=',k
1103 WRITE(*,*) ' ',valkj,val(kj), (theta/beta)**2,delta,theta
1104 END IF
1105 END IF
1106 END IF
1107 END DO ! J
1108
1109 END DO ! K
1110
1111 DO k=1,n
1112 IF(val(ilptr(k)) /= 0.0_mpd) val(ilptr(k))=1.0_mpd/val(ilptr(k))
1113 END DO
1114
1115 RETURN
1116END SUBROUTINE vabdec
1117
1123
1124SUBROUTINE vabmmm(n,val,ilptr)
1125 USE mpdef
1126
1127 IMPLICIT NONE
1128 INTEGER(mpi) :: k
1129 INTEGER(mpi) :: kp
1130 INTEGER(mpi) :: kr
1131 INTEGER(mpi) :: ks
1132 INTEGER(mpi), INTENT(IN) :: n
1133 REAL(mpd), INTENT(IN OUT) :: val(*)
1134 INTEGER(mpi), INTENT(IN) :: ilptr(n)
1135 kr=1
1136 ks=1
1137 kp=1
1138 DO k=1,n
1139 IF(val(ilptr(k)) > val(ilptr(ks))) ks=k
1140 IF(val(ilptr(k)) < val(ilptr(kr))) kr=k
1141 IF(val(ilptr(k)) > 0.0.AND.val(ilptr(k)) < val(ilptr(kp))) kp=k
1142 END DO
1143 WRITE(*,*) ' Index value ',ks,val(ilptr(ks))
1144 WRITE(*,*) ' Index value ',kp,val(ilptr(kp))
1145 WRITE(*,*) ' Index value ',kr,val(ilptr(kr))
1146
1147 RETURN
1148END SUBROUTINE vabmmm
1149
1160
1161SUBROUTINE vabslv(n,val,ilptr,x)
1162 USE mpdef
1163
1164 IMPLICIT NONE
1165 INTEGER(mpi) :: j
1166 INTEGER(mpi) :: k
1167 INTEGER(mpi) :: mk
1168
1169 INTEGER(mpi), INTENT(IN) :: n
1170 REAL(mpd), INTENT(IN OUT) :: val(*)
1171 INTEGER(mpi), INTENT(IN) :: ilptr(n)
1172 REAL(mpd), INTENT(IN OUT) :: x(n)
1173 ! ...
1174 DO k=2,n ! forward loop
1175 mk=k-ilptr(k)+ilptr(k-1)+1
1176 DO j=mk,k-1
1177 x(k)=x(k)-val(ilptr(k)-k+j)*x(j) ! X_k := X_k - L_kj B_j
1178 END DO
1179 END DO ! K
1180
1181 DO k=1,n ! divide by diagonal elements
1182 x(k)=x(k)*val(ilptr(k)) ! X_k := X_k*D_kk
1183 END DO
1184
1185 DO k=n,2,-1 ! backward loop
1186 mk=k-ilptr(k)+ilptr(k-1)+1
1187 DO j=mk,k-1
1188 x(j)=x(j)-val(ilptr(k)-k+j)*x(k) ! X_j := X_j - L_kj X_k
1189 END DO
1190 END DO ! K
1191END SUBROUTINE vabslv
1192
1193! matrix/vector products -------------------------------------------
1194
1201
1202REAL(mpd) FUNCTION dbdot(n,dx,dy)
1203 USE mpdef
1204
1205 IMPLICIT NONE
1206 INTEGER(mpi) :: i
1207 REAL(mpd) :: dtemp
1208
1209 INTEGER(mpi), INTENT(IN) :: n
1210 REAL(mpd), INTENT(IN) :: dx(n)
1211 REAL(mpd), INTENT(IN) :: dy(n)
1212 ! ...
1213 dtemp=0.0_mpd
1214 DO i = 1,mod(n,5)
1215 dtemp=dtemp+dx(i)*dy(i)
1216 END DO
1217 DO i =mod(n,5)+1,n,5
1218 dtemp=dtemp+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2) &
1219 +dx(i+3)*dy(i+3)+dx(i+4)*dy(i+4)
1220 END DO
1221 dbdot=dtemp
1222END FUNCTION dbdot
1223
1232
1233SUBROUTINE dbaxpy(n,da,dx,dy)
1234 USE mpdef
1235
1236 IMPLICIT NONE
1237 INTEGER(mpi) :: i
1238
1239 INTEGER(mpi), INTENT(IN) :: n
1240 REAL(mpd), INTENT(IN) :: dx(n)
1241 REAL(mpd), INTENT(IN OUT) :: dy(n)
1242 REAL(mpd), INTENT(IN) :: da
1243 ! ...
1244 DO i=1,mod(n,4)
1245 dy(i)=dy(i)+da*dx(i)
1246 END DO
1247 DO i=mod(n,4)+1,n,4
1248 dy(i )=dy(i )+da*dx(i )
1249 dy(i+1)=dy(i+1)+da*dx(i+1)
1250 dy(i+2)=dy(i+2)+da*dx(i+2)
1251 dy(i+3)=dy(i+3)+da*dx(i+3)
1252 END DO
1253END SUBROUTINE dbaxpy
1254
1263
1264SUBROUTINE dbsvx(v,a,b,n) !
1265 USE mpdef
1266
1267 IMPLICIT NONE
1268 INTEGER(mpi) :: i
1269 INTEGER(mpi) :: ij
1270 INTEGER(mpi) :: ijs
1271 INTEGER(mpi) :: j
1272
1273 ! B := V * A
1274 ! N N*N N
1275
1276 INTEGER(mpi), INTENT(IN) :: n
1277 REAL(mpd), INTENT(IN) :: v((n*n+n)/2)
1278 REAL(mpd), INTENT(IN) :: a(n)
1279 REAL(mpd), INTENT(OUT) :: b(n)
1280
1281 REAL(mpd) :: dsum
1282 ijs=1
1283 DO i=1,n
1284 dsum=0.0
1285 ij=ijs
1286 DO j=1,n
1287 dsum=dsum+v(ij)*a(j)
1288 IF(j < i) THEN
1289 ij=ij+1
1290 ELSE
1291 ij=ij+j
1292 END IF
1293 END DO
1294 b(i)=dsum
1295 ijs=ijs+i
1296 END DO
1297END SUBROUTINE dbsvx
1298
1307
1308SUBROUTINE dbsvxl(v,a,b,n) ! LARGE symm. matrix, vector
1309 USE mpdef
1310
1311 IMPLICIT NONE
1312 INTEGER(mpi) :: i
1313 INTEGER(mpi) :: j
1314
1315 ! B := V * A
1316 ! N N*N N
1317
1318 INTEGER(mpi), INTENT(IN) :: n
1319 REAL(mpd), INTENT(IN) :: v((INT(n,mpl)*INT(n,mpl)+INT(n,mpl))/2)
1320 REAL(mpd), INTENT(IN) :: a(n)
1321 REAL(mpd), INTENT(OUT) :: b(n)
1322
1323 REAL(mpd) :: dsum
1324 INTEGER(mpl) :: ij
1325 INTEGER(mpl) :: ijs
1326 ijs=1
1327 DO i=1,n
1328 dsum=0.0
1329 ij=ijs
1330 DO j=1,n
1331 dsum=dsum+v(ij)*a(j)
1332 IF(j < i) THEN
1333 ij=ij+1
1334 ELSE
1335 ij=ij+int(j,mpl)
1336 END IF
1337 END DO
1338 b(i)=dsum
1339 ijs=ijs+int(i,mpl)
1340 END DO
1341END SUBROUTINE dbsvxl
1342
1350
1351SUBROUTINE dbgax(a,x,y,m,n)
1352 USE mpdef
1353
1354 IMPLICIT NONE
1355 INTEGER(mpi) :: i
1356 INTEGER(mpi) :: ij
1357 INTEGER(mpi) :: j
1358
1359 INTEGER(mpi), INTENT(IN) :: m
1360 INTEGER(mpi), INTENT(IN) :: n
1361 REAL(mpd), INTENT(IN) :: a(n*m)
1362 REAL(mpd), INTENT(IN) :: x(n)
1363 REAL(mpd), INTENT(OUT) :: y(m)
1364
1365 ! ...
1366 ij=0
1367 DO i=1,m
1368 y(i)=0.0_mpd
1369 DO j=1,n
1370 ij=ij+1
1371 y(i)=y(i)+a(ij)*x(j)
1372 END DO
1373 END DO
1374END SUBROUTINE dbgax
1375
1389SUBROUTINE dbavat(v,a,w,n,m,iopt)
1390 USE mpdef
1391
1392 IMPLICIT NONE
1393 INTEGER(mpi) :: i
1394 INTEGER(mpi) :: ij
1395 INTEGER(mpi) :: ijs
1396 INTEGER(mpi) :: il
1397 INTEGER(mpi) :: j
1398 INTEGER(mpi) :: jk
1399 INTEGER(mpi) :: k
1400 INTEGER(mpi) :: l
1401 INTEGER(mpi) :: lk
1402 INTEGER(mpi) :: lkl
1403
1404 INTEGER(mpi), INTENT(IN) :: n
1405 INTEGER(mpi), INTENT(IN) :: m
1406 REAL(mpd), INTENT(IN) :: v((n*n+n)/2)
1407 REAL(mpd), INTENT(IN) :: a(n*m)
1408 REAL(mpd), INTENT(INOUT) :: w((m*m+m)/2)
1409
1410 INTEGER(mpi), INTENT(IN) :: iopt
1411
1412 REAL(mpd) :: cik
1413 ! ...
1414 IF (iopt == 0) THEN
1415 DO i=1,(m*m+m)/2
1416 w(i)=0.0_mpd ! reset output matrix
1417 END DO
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
1447END SUBROUTINE dbavat
1448
1470SUBROUTINE dbavats(v,a,is,w,n,m,iopt,sc)
1471 USE mpdef
1472
1473 IMPLICIT NONE
1474 INTEGER(mpi) :: i
1475 INTEGER(mpi) :: ic
1476 INTEGER(mpi) :: ij
1477 INTEGER(mpi) :: ijs
1478 INTEGER(mpi) :: in
1479 INTEGER(mpi) :: ir
1480 INTEGER(mpi) :: j
1481 INTEGER(mpi) :: k
1482 INTEGER(mpi) :: l
1483 INTEGER(mpi) :: lk
1484
1485 INTEGER(mpi), INTENT(IN) :: n
1486 INTEGER(mpi), INTENT(IN) :: m
1487 REAL(mpd), INTENT(IN) :: v((n*n+n)/2)
1488 REAL(mpd), INTENT(IN) :: a(n*m)
1489 REAL(mpd), INTENT(INOUT) :: w((m*m+m)/2)
1490 INTEGER(mpi), INTENT(IN) :: is(2*n*m+n+m+1)
1491 INTEGER(mpi), INTENT(IN) :: iopt
1492 INTEGER(mpi), INTENT(OUT) :: sc(n)
1493
1494 REAL(mpd) :: cik
1495 ! ...
1496 IF (iopt == 0) THEN
1497 DO i=1,(m*m+m)/2
1498 w(i)=0.0_mpd ! reset output matrix
1499 END DO
1500 END IF
1501
1502 ! offsets in V
1503 sc(1)=0
1504 DO k=2,n
1505 sc(k)=sc(k-1)+k-1
1506 END DO
1507
1508 ijs=0
1509 DO i=1,m
1510 ijs=ijs+i-1
1511 DO k=1,n
1512 cik=0.0_mpd
1513 DO l=is(i)+1,is(i+1),2
1514 ic=is(l)
1515 in=is(l+1)
1516 lk=sc(max(k,ic))+min(k,ic)
1517 cik=cik+a(in)*v(lk)
1518 END DO
1519 DO j=is(m+k)+1,is(m+k+1),2
1520 ir=is(j)
1521 in=is(j+1)
1522 IF (ir > i) EXIT
1523 ij=ijs+ir
1524 w(ij)=w(ij)+cik*a(in)
1525 END DO
1526 END DO
1527 END DO
1528END SUBROUTINE dbavats
1529
1539
1540SUBROUTINE dbmprv(lun,x,v,n)
1541 USE mpdef
1542
1543 IMPLICIT NONE
1544 INTEGER(mpi) :: i
1545 INTEGER(mpi) :: ii
1546 INTEGER(mpi) :: ij
1547 INTEGER(mpi) :: j
1548 INTEGER(mpi) :: jj
1549 INTEGER(mpi) :: l
1550 INTEGER(mpi) :: m
1551 INTEGER(mpi) :: mc(15)
1552 REAL(mps) :: pd
1553 REAL(mps) :: rho
1554 REAL(mps) :: err
1555
1556 INTEGER(mpi), INTENT(IN) :: lun
1557 INTEGER(mpi), INTENT(IN) :: n
1558 REAL(mpd), INTENT(IN) :: x(n)
1559 REAL(mpd), INTENT(IN) :: v((n*n+n)/2)
1560
1561 WRITE(lun,103)
1562 WRITE(lun,101)
1563 ii=0
1564 DO i=1,n
1565 ij=ii
1566 ii=ii+i
1567 err=0.0
1568 IF(v(ii) > 0.0) err=sqrt(real(v(ii),mps))
1569 l=0
1570 jj=0
1571 DO j=1,i
1572 jj=jj+j
1573 ij=ij+1
1574 rho=0.0
1575 pd=real(v(ii)*v(jj),mps)
1576 IF(pd > 0.0) rho=real(v(ij),mps)/sqrt(pd)
1577 l=l+1
1578 mc(l)=nint(100.0*abs(rho),mpi)
1579 IF(rho < 0.0) mc(l)=-mc(l)
1580 IF(j == i.OR.l == 15) THEN
1581 IF(j <= 15) THEN
1582 IF(j == i) THEN
1583 WRITE(lun,102) i,x(i),err,(mc(m),m=1,l-1)
1584 ELSE
1585 WRITE(lun,102) i,x(i),err,(mc(m),m=1,l)
1586 END IF
1587 ELSE
1588 IF(j == i) THEN
1589 WRITE(lun,103) (mc(m),m=1,l-1)
1590 ELSE
1591 WRITE(lun,103) (mc(m),m=1,l)
1592 END IF
1593 l=0
1594 END IF
1595 END IF
1596 END DO
1597 END DO
1598 WRITE(lun,104)
1599 ! 100 RETURN
1600 RETURN
1601101 FORMAT(9x,'Param',7x,'error',7x,'correlation coefficients'/)
1602102 FORMAT(1x,i5,2g12.4,1x,15i5)
1603103 FORMAT(31x,15i5)
1604104 FORMAT(33x,'(correlation coefficients in percent)')
1605END SUBROUTINE dbmprv
1606
1614
1615SUBROUTINE dbprv(lun,v,n)
1616 USE mpdef
1617
1618 IMPLICIT NONE
1619 INTEGER(mpi), PARAMETER :: istp=6
1620 INTEGER(mpi) :: i
1621 INTEGER(mpi) :: ip
1622 INTEGER(mpi) :: ipe
1623 INTEGER(mpi) :: ipn
1624 INTEGER(mpi) :: ips
1625 INTEGER(mpi) :: k
1626
1627 INTEGER(mpi), INTENT(IN) :: lun
1628 INTEGER(mpi), INTENT(IN) :: n
1629 REAL(mpd), INTENT(IN) :: v((n*n+n)/2)
1630
1631 WRITE(lun,101)
1632
1633 DO i=1,n
1634 ips=(i*i-i)/2
1635 ipe=ips+i
1636 ip =ips
1637100 CONTINUE
1638 ipn=ip+istp
1639 WRITE(lun,102) i, ip+1-ips, (v(k),k=ip+1,min(ipn,ipe))
1640 IF (ipn < ipe) THEN
1641 ip=ipn
1642 GO TO 100
1643 END IF
1644END DO
1645RETURN
1646101 FORMAT(1x,'--- DBPRV -----------------------------------')
1647102 FORMAT(1x,2i3,6g12.4)
1648END SUBROUTINE dbprv
1649
1650! sort -------------------------------------------------------------
1651
1658
1659SUBROUTINE heapf(a,n)
1660 USE mpdef
1661
1662 IMPLICIT NONE
1663 !
1664 INTEGER(mpi) :: i
1665 INTEGER(mpi) :: j
1666 INTEGER(mpi) :: l
1667 INTEGER(mpi) :: r
1668 REAL(mps) :: at ! pivot key value
1669
1670 INTEGER(mpi), INTENT(IN) :: n
1671 REAL(mps), INTENT(IN OUT) :: a(n)
1672 ! ...
1673 IF(n <= 1) RETURN
1674 l=n/2+1
1675 r=n
167610 IF(l > 1) THEN
1677 l=l-1
1678 at =a(l)
1679 ELSE
1680 at =a(r)
1681 a(r)=a(1)
1682 r=r-1
1683 IF(r == 1) THEN
1684 a(1)=at
1685 RETURN
1686 END IF
1687 END IF
1688 i=l
1689 j=l+l
169020 IF(j <= r) THEN
1691 IF(j < r) THEN
1692 IF(a(j) < a(j+1)) j=j+1
1693 END IF
1694 IF(at < a(j)) THEN
1695 a(i)=a(j)
1696 i=j
1697 j=j+j
1698 ELSE
1699 j=r+1
1700 END IF
1701 GO TO 20
1702 END IF
1703 a(i)=at
1704 GO TO 10
1705END SUBROUTINE heapf
1706
1713
1714SUBROUTINE sort1k(a,n)
1715 USE mpdef
1716
1717 IMPLICIT NONE
1718 INTEGER(mpi) :: nlev ! stack size
1719 parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1720 INTEGER(mpi) :: i
1721 INTEGER(mpi) :: j
1722 INTEGER(mpi) :: l
1723 INTEGER(mpi) :: r
1724 INTEGER(mpi) :: lev
1725 INTEGER(mpi) :: lr(nlev)
1726 INTEGER(mpi) :: lrh
1727 INTEGER(mpi) :: maxlev
1728 INTEGER(mpi) :: a1 ! pivot key
1729 INTEGER(mpi) :: at ! pivot key
1730
1731 INTEGER(mpi), INTENT(IN) :: n
1732 INTEGER(mpi), INTENT(IN OUT) :: a(n)
1733 ! ...
1734 IF (n <= 0) RETURN
1735 maxlev=0
1736 lev=0
1737 l=1
1738 r=n
173910 IF(r-l == 1) THEN ! sort two elements L and R
1740 IF(a(l) > a(r)) THEN
1741 at=a(l) ! exchange L <-> R
1742 a(l)=a(r)
1743 a(r)=at
1744 END IF
1745 r=l
1746 END IF
1747 IF(r == l) THEN
1748 IF(lev <= 0) THEN
1749 ! WRITE(*,*) 'SORT1K (quicksort): maxlevel used/available =',
1750 ! + MAXLEV,'/64'
1751 RETURN
1752 END IF
1753 lev=lev-2
1754 l=lr(lev+1)
1755 r=lr(lev+2)
1756 ELSE
1757 ! LRH=(L+R)/2
1758 lrh=(l/2)+(r/2) ! avoid bit overflow
1759 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1760 a1=a(lrh) ! middle
1761 i=l-1 ! find limits [J,I] with [L,R]
1762 j=r+1
176320 i=i+1
1764 IF(a(i) < a1) GO TO 20
176530 j=j-1
1766 IF(a(j) > a1) GO TO 30
1767 IF(i <= j) THEN
1768 at=a(i) ! exchange I <-> J
1769 a(i)=a(j)
1770 a(j)=at
1771 GO TO 20
1772 END IF
1773 IF(lev+2 > nlev) THEN
1774 CALL peend(33,'Aborted, stack overflow in quicksort')
1775 stop 'SORT1K (quicksort): stack overflow'
1776 END IF
1777 IF(r-i < j-l) THEN
1778 lr(lev+1)=l
1779 lr(lev+2)=j
1780 l=i
1781 ELSE
1782 lr(lev+1)=i
1783 lr(lev+2)=r
1784 r=j
1785 END IF
1786 lev=lev+2
1787 maxlev=max(maxlev,lev)
1788 END IF
1789 GO TO 10
1790END SUBROUTINE sort1k
1791
1798
1799SUBROUTINE sort2k(a,n)
1800 USE mpdef
1801
1802 IMPLICIT NONE
1803 INTEGER(mpi) :: nlev ! stack size
1804 parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1805 INTEGER(mpi) :: i
1806 INTEGER(mpi) ::j
1807 INTEGER(mpi) ::l
1808 INTEGER(mpi) ::r
1809 INTEGER(mpi) ::lev
1810 INTEGER(mpi) ::lr(nlev)
1811 INTEGER(mpi) ::lrh
1812 INTEGER(mpi) ::maxlev
1813 INTEGER(mpi) ::a1 ! pivot key
1814 INTEGER(mpi) ::a2 ! pivot key
1815 INTEGER(mpi) ::at ! pivot key
1816
1817 INTEGER(mpi), INTENT(IN OUT) :: a(2,*)
1818 INTEGER(mpi), INTENT(IN) :: n
1819 ! ...
1820 maxlev=0
1821 lev=0
1822 l=1
1823 r=n
182410 IF(r-l == 1) THEN ! sort two elements L and R
1825 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
1826 at=a(1,l) ! exchange L <-> R
1827 a(1,l)=a(1,r)
1828 a(1,r)=at
1829 at=a(2,l)
1830 a(2,l)=a(2,r)
1831 a(2,r)=at
1832 END IF
1833 r=l
1834 END IF
1835 IF(r == l) THEN
1836 IF(lev <= 0) THEN
1837 WRITE(*,*) 'SORT2K (quicksort): maxlevel used/available =', maxlev,'/64'
1838 RETURN
1839 END IF
1840 lev=lev-2
1841 l=lr(lev+1)
1842 r=lr(lev+2)
1843 ELSE
1844 ! LRH=(L+R)/2
1845 lrh=(l/2)+(r/2) ! avoid bit overflow
1846 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1847 a1=a(1,lrh) ! middle
1848 a2=a(2,lrh)
1849 i=l-1 ! find limits [J,I] with [L,R]
1850 j=r+1
185120 i=i+1
1852 IF(a(1,i) < a1) GO TO 20
1853 IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
185430 j=j-1
1855 IF(a(1,j) > a1) GO TO 30
1856 IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
1857 IF(i <= j) THEN
1858 at=a(1,i) ! exchange I <-> J
1859 a(1,i)=a(1,j)
1860 a(1,j)=at
1861 at=a(2,i)
1862 a(2,i)=a(2,j)
1863 a(2,j)=at
1864 GO TO 20
1865 END IF
1866 IF(lev+2 > nlev) THEN
1867 CALL peend(33,'Aborted, stack overflow in quicksort')
1868 stop 'SORT2K (quicksort): stack overflow'
1869 END IF
1870 IF(r-i < j-l) THEN
1871 lr(lev+1)=l
1872 lr(lev+2)=j
1873 l=i
1874 ELSE
1875 lr(lev+1)=i
1876 lr(lev+2)=r
1877 r=j
1878 END IF
1879 lev=lev+2
1880 maxlev=max(maxlev,lev)
1881 END IF
1882 GO TO 10
1883END SUBROUTINE sort2k
1884
1891
1892SUBROUTINE sort2i(a,n)
1893 USE mpdef
1894
1895 IMPLICIT NONE
1896 INTEGER(mpi) :: nlev ! stack size
1897 parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1898 INTEGER(mpi) :: i
1899 INTEGER(mpi) ::j
1900 INTEGER(mpi) ::l
1901 INTEGER(mpi) ::r
1902 INTEGER(mpi) ::lev
1903 INTEGER(mpi) ::lr(nlev)
1904 INTEGER(mpi) ::lrh
1905 INTEGER(mpi) ::maxlev
1906 INTEGER(mpi) ::a1 ! pivot key
1907 INTEGER(mpi) ::a2 ! pivot key
1908 INTEGER(mpi) ::at(3)
1909
1910 INTEGER(mpi), INTENT(IN OUT) :: a(3,*)
1911 INTEGER(mpi), INTENT(IN) :: n
1912 ! ...
1913 maxlev=0
1914 lev=0
1915 l=1
1916 r=n
191710 IF(r-l == 1) THEN ! sort two elements L and R
1918 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
1919 at=a(:,l) ! exchange L <-> R
1920 a(:,l)=a(:,r)
1921 a(:,r)=at
1922 END IF
1923 r=l
1924 END IF
1925 IF(r == l) THEN
1926 IF(lev <= 0) THEN
1927 WRITE(*,*) 'SORT2I (quicksort): maxlevel used/available =', maxlev,'/64'
1928 RETURN
1929 END IF
1930 lev=lev-2
1931 l=lr(lev+1)
1932 r=lr(lev+2)
1933 ELSE
1934 ! LRH=(L+R)/2
1935 lrh=(l/2)+(r/2) ! avoid bit overflow
1936 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1937 a1=a(1,lrh) ! middle
1938 a2=a(2,lrh)
1939 i=l-1 ! find limits [J,I] with [L,R]
1940 j=r+1
194120 i=i+1
1942 IF(a(1,i) < a1) GO TO 20
1943 IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
194430 j=j-1
1945 IF(a(1,j) > a1) GO TO 30
1946 IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
1947 IF(i <= j) THEN
1948 IF(a(1,i) == a(1,j).AND.a(2,i) == a(2,j)) GO TO 20 ! equal -> keep order
1949 at=a(:,i) ! exchange I <-> J
1950 a(:,i)=a(:,j)
1951 a(:,j)=at
1952 GO TO 20
1953 END IF
1954 IF(lev+2 > nlev) THEN
1955 CALL peend(33,'Aborted, stack overflow in quicksort')
1956 stop 'SORT2I (quicksort): stack overflow'
1957 END IF
1958 IF(r-i < j-l) THEN
1959 lr(lev+1)=l
1960 lr(lev+2)=j
1961 l=i
1962 ELSE
1963 lr(lev+1)=i
1964 lr(lev+2)=r
1965 r=j
1966 END IF
1967 lev=lev+2
1968 maxlev=max(maxlev,lev)
1969 END IF
1970 GO TO 10
1971END SUBROUTINE sort2i
1972
1980
1981SUBROUTINE sort22l(a,b,n)
1982 USE mpdef
1983
1984 IMPLICIT NONE
1985 INTEGER(mpi) :: nlev ! stack size
1986 parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1987 INTEGER(mpi) :: i
1988 INTEGER(mpi) ::j
1989 INTEGER(mpi) ::l
1990 INTEGER(mpi) ::r
1991 INTEGER(mpi) ::lev
1992 INTEGER(mpi) ::lr(nlev)
1993 INTEGER(mpi) ::lrh
1994 INTEGER(mpi) ::maxlev
1995 INTEGER(mpi) ::a1 ! pivot key
1996 INTEGER(mpi) ::a2 ! pivot key
1997 INTEGER(mpi) ::at(4)
1998 INTEGER(mpl) ::bt
1999
2000 INTEGER(mpi), INTENT(IN OUT) :: a(4,*)
2001 INTEGER(mpl), INTENT(IN OUT) :: b(*)
2002 INTEGER(mpi), INTENT(IN) :: n
2003 ! ...
2004 maxlev=0
2005 lev=0
2006 l=1
2007 r=n
2008 IF (n<=1) RETURN
200910 IF(r-l == 1) THEN ! sort two elements L and R
2010 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
2011 at=a(:,l) ! exchange L <-> R
2012 a(:,l)=a(:,r)
2013 a(:,r)=at
2014 bt=b(l)
2015 b(l)=b(r)
2016 b(r)=bt
2017 END IF
2018 r=l
2019 END IF
2020 IF(r == l) THEN
2021 IF(lev <= 0) THEN
2022 WRITE(*,*) 'SORT22l (quicksort): maxlevel used/available =', maxlev,'/64'
2023 RETURN
2024 END IF
2025 lev=lev-2
2026 l=lr(lev+1)
2027 r=lr(lev+2)
2028 ELSE
2029 ! LRH=(L+R)/2
2030 lrh=(l/2)+(r/2) ! avoid bit overflow
2031 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
2032 a1=a(1,lrh) ! middle
2033 a2=a(2,lrh)
2034 i=l-1 ! find limits [J,I] with [L,R]
2035 j=r+1
203620 i=i+1
2037 IF(a(1,i) < a1) GO TO 20
2038 IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
203930 j=j-1
2040 IF(a(1,j) > a1) GO TO 30
2041 IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
2042 IF(i <= j) THEN
2043 at=a(:,i) ! exchange I <-> J
2044 a(:,i)=a(:,j)
2045 a(:,j)=at
2046 bt=b(i)
2047 b(i)=b(j)
2048 b(j)=bt
2049 GO TO 20
2050 END IF
2051 IF(lev+2 > nlev) THEN
2052 CALL peend(33,'Aborted, stack overflow in quicksort')
2053 stop 'SORT22l (quicksort): stack overflow'
2054 END IF
2055 IF(r-i < j-l) THEN
2056 lr(lev+1)=l
2057 lr(lev+2)=j
2058 l=i
2059 ELSE
2060 lr(lev+1)=i
2061 lr(lev+2)=r
2062 r=j
2063 END IF
2064 lev=lev+2
2065 maxlev=max(maxlev,lev)
2066 END IF
2067 GO TO 10
2068END SUBROUTINE sort22l
2069
2077
2078REAL(mps) FUNCTION chindl(n,nd)
2079 USE mpdef
2080
2081 IMPLICIT NONE
2082 INTEGER(mpi) :: m
2083 INTEGER(mpi), INTENT(IN) :: n
2084 INTEGER(mpi), INTENT(IN) :: nd
2085 !
2086 REAL(mps) :: sn(3)
2087 REAL(mps) ::table(30,3)
2088 ! REAL PN(3)
2089 ! DATA PN/0.31731,0.0455002785,2.69985E-3/ ! probabilities
2090 DATA sn/0.47523,1.690140,2.782170/
2091 DATA table/ 1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, &
2092 1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, &
2093 1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, &
2094 1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040, &
2095 4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, &
2096 1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, &
2097 1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, &
2098 1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742, &
2099 9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, &
2100 2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, &
2101 2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, &
2102 1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681/
2103 SAVE sn,table
2104 ! ...
2105 IF(nd < 1) THEN
2106 chindl=0.0
2107 ELSE
2108 m=max(1,min(n,3)) ! 1, 2 or 3 sigmas
2109 IF(nd <= 30) THEN
2110 chindl=table(nd,m) ! from table
2111 ELSE ! approximation for ND > 30
2112 chindl=(sn(m)+sqrt(real(nd+nd-1,mps)))**2/real(nd+nd,mps)
2113 END IF
2114 END IF
2115END FUNCTION chindl
2116
2154
2155SUBROUTINE lltdec(n,c,india,nrkd,iopt)
2156 USE mpdef
2157
2158 IMPLICIT NONE
2159 INTEGER(mpi) :: i
2160 INTEGER(mpi) :: j
2161 INTEGER(mpi) :: k
2162 INTEGER(mpi) :: kj
2163 INTEGER(mpi) :: mj
2164 INTEGER(mpi) :: mk
2165 REAL(mpd) ::diag
2166
2167 INTEGER(mpi), INTENT(IN) :: n
2168 INTEGER(mpi), INTENT(IN) :: india(n)
2169 INTEGER(mpi), INTENT(OUT) :: nrkd
2170 INTEGER(mpi), INTENT(IN) :: iopt
2171 REAL(mpd), INTENT(IN OUT) :: c(india(n))
2172 REAL(mpd) eps
2173 ! ...
2174 eps = 16.0_mpd * epsilon(eps) ! 16 * precision(mpd)
2175
2176 ! ..
2177 nrkd=0
2178 diag=0.0_mpd
2179 IF(c(india(1)) > 0.0) THEN
2180 c(india(1))=1.0_mpd/sqrt(c(india(1))) ! square root
2181 ELSE
2182 c(india(1))=0.0_mpd
2183 nrkd=-1
2184 END IF
2185
2186 DO k=2,n
2187 mk=k-india(k)+india(k-1)+1 ! first index in row K
2188 DO j=mk,k ! loop over row K with index J
2189 mj=1
2190 IF(j > 1) mj=j-india(j)+india(j-1)+1 ! first index in row J
2191 kj=india(k)-k+j ! index kj
2192 diag=c(india(j)) ! j-th diagonal element
2193
2194 DO i=max(mj,mk),j-1
2195 ! L_kj = L_kj - L_ki *D_ii *L_ji
2196 c(kj)=c(kj) - c(india(k)-k+i)*c(india(j)-j+i)
2197 END DO ! I
2198
2199 IF(j /= k) c(kj)=c(kj)*diag
2200 END DO ! J
2201
2202 IF(c(india(k)) > eps*diag) THEN ! test for linear dependence
2203 c(india(k))=1.0_mpd/sqrt(c(india(k))) ! square root
2204 ELSE
2205 DO j=mk,k ! reset row K
2206 c(india(k)-k+j)=0.0_mpd
2207 END DO ! J
2208 IF (iopt > 0 .and. diag > 0.0) THEN ! skyline
2209 c(india(k))=1.0_mpd/sqrt(diag) ! square root
2210 ELSE
2211 IF(nrkd == 0) THEN
2212 nrkd=-k
2213 ELSE
2214 IF(nrkd < 0) nrkd=1
2215 nrkd=nrkd+1
2216 END IF
2217 END IF
2218 END IF
2219
2220 END DO ! K
2221 RETURN
2222END SUBROUTINE lltdec
2223
2224
2236
2237SUBROUTINE lltfwd(n,c,india,x)
2238 USE mpdef
2239
2240 IMPLICIT NONE
2241 INTEGER(mpi) :: j
2242 INTEGER(mpi) :: k
2243
2244 INTEGER(mpi), INTENT(IN) :: n
2245 INTEGER(mpi), INTENT(IN) :: india(n)
2246 REAL(mpd), INTENT(IN OUT) :: c(india(n))
2247 REAL(mpd), INTENT(IN OUT) :: x(n)
2248
2249 x(1)=x(1)*c(india(1))
2250 DO k=2,n ! forward loop
2251 DO j=k-india(k)+india(k-1)+1,k-1
2252 x(k)=x(k)-c(india(k)-k+j)*x(j) ! X_k := X_k - L_kj * B_j
2253 END DO ! J
2254 x(k)=x(k)*c(india(k))
2255 END DO ! K
2256
2257 RETURN
2258END SUBROUTINE lltfwd
2259
2274
2275SUBROUTINE lltfwds(n,c,india,x,i0,ns)
2276 USE mpdef
2277
2278 IMPLICIT NONE
2279 INTEGER(mpi) :: j
2280 INTEGER(mpi) :: k
2281 INTEGER(mpi) :: l
2282
2283 INTEGER(mpi), INTENT(IN) :: n
2284 INTEGER(mpi), INTENT(IN) :: india(n)
2285 REAL(mpd), INTENT(IN OUT) :: c(india(n))
2286 REAL(mpd), INTENT(IN OUT) :: x(n)
2287 INTEGER(mpi), INTENT(IN) :: i0
2288 INTEGER(mpi), INTENT(IN) :: ns
2289
2290 DO k=1,min(ns,n-i0) ! forward loop
2291 l=k+i0
2292 IF (l>1) THEN
2293 DO j=max(1,k-india(l)+india(l-1)+1),k-1
2294 x(k)=x(k)-c(india(l)-k+j)*x(j) ! X_k := X_k - L_kj * B_j
2295 END DO ! J
2296 END IF
2297 x(k)=x(k)*c(india(l))
2298 END DO ! K
2299
2300 RETURN
2301END SUBROUTINE lltfwds
2302
2314
2315SUBROUTINE lltbwd(n,c,india,x)
2316 USE mpdef
2317
2318 IMPLICIT NONE
2319 INTEGER(mpi) :: j
2320 INTEGER(mpi) :: k
2321
2322 INTEGER(mpi), INTENT(IN) :: n
2323 INTEGER(mpi), INTENT(IN) :: india(n)
2324 REAL(mpd), INTENT(IN OUT) :: c(india(n))
2325 REAL(mpd), INTENT(IN OUT) :: x(n)
2326
2327 DO k=n,2,-1 ! backward loop
2328 x(k)=x(k)*c(india(k))
2329 DO j=k-india(k)+india(k-1)+1,k-1
2330 x(j)=x(j)-c(india(k)-k+j)*x(k) ! X_j := X_j - L_kj * X_k
2331 END DO ! J
2332 END DO ! K
2333 x(1)=x(1)*c(india(1))
2334
2335 RETURN
2336END SUBROUTINE lltbwd
2337
2355SUBROUTINE equdec(n,m,ls,c,india,nrkd,nrkd2)
2356 USE mpdef
2357
2358 IMPLICIT NONE
2359 INTEGER(mpi) :: i
2360 INTEGER(mpi) :: j
2361 INTEGER(mpl) :: jk
2362 INTEGER(mpi) :: k
2363 INTEGER(mpl) :: nn
2364
2365 INTEGER(mpi), INTENT(IN) :: n
2366 INTEGER(mpi), INTENT(IN) :: m
2367 INTEGER(mpi), INTENT(IN) :: ls
2368 INTEGER(mpi), INTENT(IN OUT) :: india(n+m)
2369 INTEGER(mpi), INTENT(OUT) :: nrkd
2370 INTEGER(mpi), INTENT(OUT) :: nrkd2
2371 REAL(mpd), INTENT(IN OUT) :: c(india(n)+n*m+(m*m+m)/2)
2372
2373 ! ...
2374
2375 nrkd=0
2376 nrkd2=0
2377 nn=n
2378
2379 CALL lltdec(n,c,india,nrkd,ls) ! decomposition G G^T
2380
2381 IF (m>0) THEN
2382 DO i=1,m
2383 CALL lltfwd(n,c,india,c(india(n)+int(i-1,mpl)*nn+1)) ! forward solution K
2384 END DO
2385
2386 jk=india(n)+nn*int(m,mpl)
2387 DO j=1,m
2388 DO k=1,j
2389 jk=jk+1
2390 c(jk)=0.0_mpd ! product K K^T
2391 DO i=1,n
2392 c(jk)=c(jk)+c(india(n)+int(j-1,mpl)*nn+i)*c(india(n)+int(k-1,mpl)*nn+i)
2393 END DO
2394 END DO
2395 END DO
2396
2397 india(n+1)=1
2398 DO i=2,m
2399 india(n+i)=india(n+i-1)+min(i,m) ! pointer for K K^T
2400 END DO
2401 CALL lltdec(m,c(india(n)+nn*int(m,mpl)+1),india(n+1),nrkd2,0) ! decomp. H H^T
2402 ENDIF
2403
2404 RETURN
2405END SUBROUTINE equdec
2406
2422SUBROUTINE equslv(n,m,c,india,x) ! solution vector
2423 USE mpdef
2424
2425 IMPLICIT NONE
2426 INTEGER(mpi) :: i
2427 INTEGER(mpi) :: j
2428 INTEGER(mpl) :: nn
2429
2430 INTEGER(mpi), INTENT(IN) :: n
2431 INTEGER(mpi), INTENT(IN) :: m
2432 INTEGER(mpi), INTENT(IN) :: india(n+m)
2433 REAL(mpd), INTENT(IN OUT) :: c(india(n)+n*m+(m*m+m)/2)
2434 REAL(mpd), INTENT(IN OUT) :: x(n+m)
2435
2436 CALL lltfwd(n,c,india,x) ! result is u
2437
2438 IF (m>0) THEN
2439 nn=n
2440 DO i=1,m
2441 DO j=1,n
2442 x(n+i)=x(n+i)-x(j)*c(india(n)+int(i-1,mpl)*nn+j) ! g - K u
2443 END DO
2444 END DO
2445 CALL lltfwd(m,c(india(n)+nn*int(m,mpl)+1),india(n+1),x(n+1)) ! result is v
2446
2447
2448 CALL lltbwd(m,c(india(n)+nn*int(m,mpl)+1),india(n+1),x(n+1)) ! result is -y
2449 DO i=1,m
2450 x(n+i)=-x(n+i) ! result is +y
2451 END DO
2452
2453 DO i=1,n
2454 DO j=1,m
2455 x(i)=x(i)-x(n+j)*c(india(n)+int(j-1,mpl)*nn+i) ! u - K^T y
2456 END DO
2457 END DO
2458 ENDIF
2459
2460 CALL lltbwd(n,c,india,x) ! result is x
2461
2462 RETURN
2463END SUBROUTINE equslv
2464
2486SUBROUTINE equdecs(n,m,b,nm,ls,c,india,l,nrkd,nrkd2)
2487 USE mpdef
2488
2489 IMPLICIT NONE
2490 INTEGER(mpi) :: i
2491 INTEGER(mpi) :: i0
2492 INTEGER(mpi) :: ib
2493 INTEGER(mpi) :: in
2494 INTEGER(mpi) :: j
2495 INTEGER(mpi) :: j0
2496 INTEGER(mpi) :: jb
2497 INTEGER(mpi) :: jk
2498 INTEGER(mpi) :: jn
2499 INTEGER(mpi) :: k
2500 INTEGER(mpi) :: nc
2501 INTEGER(mpi) :: p1
2502 INTEGER(mpi) :: p2
2503 INTEGER(mpi) :: q1
2504 INTEGER(mpi) :: q2
2505
2506 INTEGER(mpl) :: ij
2507 INTEGER(mpl) :: jj
2508 INTEGER(mpl) :: kk
2509
2510 INTEGER(mpi), INTENT(IN) :: n
2511 INTEGER(mpi), INTENT(IN) :: m
2512 INTEGER(mpi), INTENT(IN) :: b
2513 INTEGER(mpl), INTENT(IN) :: nm
2514 INTEGER(mpi), INTENT(IN) :: ls
2515 INTEGER(mpi), INTENT(IN OUT) :: india(n+m)
2516 INTEGER(mpi), INTENT(IN) :: l(4*b)
2517 INTEGER(mpi), INTENT(OUT) :: nrkd
2518 INTEGER(mpi), INTENT(OUT) :: nrkd2
2519 REAL(mpd), INTENT(IN OUT) :: c(nm+india(n)+(m*m+m)/2)
2520
2521 ! ...
2522
2523 nrkd=0
2524 nrkd2=0
2525
2526 CALL lltdec(n,c,india,nrkd,ls) ! decomposition G G^T
2527
2528 IF (m>0) THEN
2529 ij=india(n)+(m*m+m)/2 ! offset of block part
2530 DO ib=1,b ! loop over blocks
2531 p1=l(ib*4-3) ! first constraint in block
2532 p2=l(ib*4-2) ! last constraint in block
2533 i0=l(ib*4-1) ! parameter offset in block
2534 in=l(ib*4) ! number of columns in block
2535 DO j=p1,p2
2536 CALL lltfwds(n,c,india,c(ij+1),i0,in) ! forward solution K
2537 ij=ij+in
2538 END DO
2539 END DO
2540
2541 ij=india(n)+(m*m+m)/2 ! offset of block part
2542 DO ib=1,b ! loop over blocks
2543 p1=l(ib*4-3) ! first constraint in block
2544 p2=l(ib*4-2) ! last constraint in block
2545 i0=l(ib*4-1) ! parameter offset in block
2546 in=l(ib*4) ! number of columns in block
2547 DO i=1,in ! loop over columns in block
2548 ij=ij+1
2549 DO j=p1,p2
2550 jk=india(n)+(j*j-j)/2+p1
2551 DO k=p1,j
2552 c(jk)=c(jk)+c(ij+(j-p1)*in)*c(ij+(k-p1)*in) ! product K K^T
2553 jk=jk+1
2554 END DO
2555 END DO
2556 END DO
2557 jj=ij+(p2-p1)*in ! offset of next block
2558 ! loop over potentially overlapping blocks
2559 DO jb=ib+1,b
2560 q1=l(jb*4-3) ! first constraint in block
2561 q2=l(jb*4-2) ! last constraint in block
2562 j0=l(jb*4-1) ! parameter offset in block
2563 jn=l(jb*4) ! number of columns in block
2564 nc=min(in+i0-j0,jn) ! overlapping columns
2565 IF (nc < 1) EXIT
2566 !print *, ' EQUDECS overlap ', ib,jb,nc,i0,in,j0,jn
2567 kk=ij+j0-i0-in
2568 DO i=1,nc ! loop over parameters in block
2569 kk=kk+1
2570 DO j=q1,q2
2571 jk=india(n)+(j*j-j)/2+p1
2572 DO k=p1,p2
2573 c(jk)=c(jk)+c(jj+(j-q1)*jn+i)*c(kk+(k-p1)*in) ! product K K^T
2574 jk=jk+1
2575 END DO
2576 END DO
2577 END DO
2578 jj=jj+(q2+1-q1)*jn
2579 END DO
2580 ij=ij+(p2-p1)*in
2581 END DO
2582
2583 india(n+1)=1
2584 DO i=2,m
2585 india(n+i)=india(n+i-1)+min(i,m) ! pointer for K K^T
2586 END DO
2587
2588 CALL lltdec(m,c(india(n)+1),india(n+1),nrkd2,0) ! decomp. H H^T
2589 ENDIF
2590
2591 RETURN
2592END SUBROUTINE equdecs
2593
2613SUBROUTINE equslvs(n,m,b,nm,c,india,l,x) ! solution vector
2614 USE mpdef
2615
2616 IMPLICIT NONE
2617 INTEGER(mpi) :: i
2618 INTEGER(mpi) :: i0
2619 INTEGER(mpi) :: ib
2620 INTEGER(mpi) :: in
2621 INTEGER(mpi) :: j
2622 INTEGER(mpi) :: p1
2623 INTEGER(mpi) :: p2
2624
2625 INTEGER(mpl) :: ij
2626
2627 INTEGER(mpi), INTENT(IN) :: n
2628 INTEGER(mpi), INTENT(IN) :: m
2629 INTEGER(mpi), INTENT(IN) :: b
2630 INTEGER(mpl), INTENT(IN) :: nm
2631 INTEGER(mpi), INTENT(IN) :: india(n+m)
2632 INTEGER(mpi), INTENT(IN) :: l(4*b)
2633 REAL(mpd), INTENT(IN OUT) :: c(nm+india(n)+(m*m+m)/2)
2634 REAL(mpd), INTENT(IN OUT) :: x(n+m)
2635
2636 CALL lltfwd(n,c,india,x) ! result is u
2637
2638 IF (m>0) THEN
2639 ij=india(n)+(m*m+m)/2 ! offset of block part
2640 DO ib=1,b ! loop over blocks
2641 p1=l(ib*4-3) ! first constraint in block
2642 p2=l(ib*4-2) ! last constraint in block
2643 i0=l(ib*4-1) ! parameter offset in block
2644 in=l(ib*4) ! number of columns in block
2645 DO j=p1,p2 ! loop over constraints in block
2646 DO i=1,in ! loop over parameters in block
2647 ij=ij+1
2648 x(n+j)=x(n+j)-c(ij)*x(i+i0) ! g - K u
2649 END DO
2650 END DO
2651 END DO
2652
2653 CALL lltfwd(m,c(india(n)+1),india(n+1),x(n+1)) ! result is v
2654
2655 CALL lltbwd(m,c(india(n)+1),india(n+1),x(n+1)) ! result is -y
2656 DO i=1,m
2657 x(n+i)=-x(n+i) ! result is +y
2658 END DO
2659
2660 ij=india(n)+(m*m+m)/2 ! offset of block part
2661 DO ib=1,b ! loop over blocks
2662 p1=l(ib*4-3) ! first constraint in block
2663 p2=l(ib*4-2) ! last constraint in block
2664 i0=l(ib*4-1) ! parameter offset in block
2665 in=l(ib*4) ! number of columns in block
2666 DO j=p1,p2 ! loop over constraints in block
2667 DO i=1,in ! loop over parameters in block
2668 ij=ij+1
2669 x(i+i0)=x(i+i0)-c(ij)*x(n+j) ! u - K^T y
2670 END DO
2671 END DO
2672 END DO
2673 ENDIF
2674
2675 CALL lltbwd(n,c,india,x) ! result is x
2676
2677 RETURN
2678END SUBROUTINE equslvs
2679
2709
2710SUBROUTINE precon(p,n,c,cu,a,s,nrkd)
2711 USE mpdef
2712
2713 IMPLICIT NONE
2714 INTEGER(mpi) :: i
2715 INTEGER(mpi) :: ii
2716 INTEGER(mpi) :: j
2717 INTEGER(mpi) :: jj
2718 INTEGER(mpi) :: jk
2719 INTEGER(mpi) :: k
2720 INTEGER(mpi) :: kk
2721
2722 INTEGER(mpi), INTENT(IN) :: p
2723 INTEGER(mpi), INTENT(IN) :: n
2724 REAL(mpd), INTENT(IN) :: c(n)
2725 REAL(mpd), INTENT(OUT) :: cu(n)
2726 REAL(mpd), INTENT(IN OUT) :: a(n,p)
2727 REAL(mpd), INTENT(OUT) :: s((p*p+p)/2)
2728 INTEGER(mpi), INTENT(OUT) :: nrkd
2729
2730 REAL(mpd) :: div
2731 REAL(mpd) :: ratio
2732
2733 nrkd=0
2734 DO i=1,(p*p+p)/2
2735 s(i)=0.0_mpd
2736 END DO
2737 DO i=1,n
2738 jk=0
2739 div=c(i) ! copy
2740 IF (div > 0.0_mpd) THEN
2741 cu(i)=1.0_mpd/sqrt(div)
2742 ELSE
2743 cu(i)=0.0_mpd
2744 nrkd=nrkd+1
2745 END IF
2746 DO j=1,p
2747 a(i,j)=a(i,j)*cu(i) ! K = A C^{-1/2}
2748 DO k=1,j
2749 jk=jk+1
2750 s(jk)=s(jk)+a(i,j)*a(i,k) ! S = symmetric matrix K K^T
2751 END DO
2752 END DO
2753 END DO
2754
2755 ii=0
2756 DO i=1,p ! S -> H D H^T (Cholesky)
2757 ii=ii+i
2758 IF(s(ii) /= 0.0_mpd) s(ii)=1.0_mpd/s(ii)
2759 jj=ii
2760 DO j=i+1,p
2761 ratio=s(i+jj)*s(ii)
2762 kk=jj
2763 DO k=j,p
2764 s(kk+j)=s(kk+j)-s(kk+i)*ratio
2765 kk=kk+k
2766 END DO ! K
2767 s(i+jj)=ratio
2768 jj=jj+j
2769 END DO ! J
2770 END DO ! I
2771 RETURN
2772END SUBROUTINE precon
2773
2783
2784SUBROUTINE presol(p,n,cu,a,s,x,y) ! solution
2785 USE mpdef
2786
2787 IMPLICIT NONE
2788 INTEGER(mpi) :: i
2789 INTEGER(mpi) :: j
2790 INTEGER(mpi) :: jj
2791 INTEGER(mpi) :: k
2792 INTEGER(mpi) :: kk
2793
2794 INTEGER(mpi), INTENT(IN) :: p
2795 INTEGER(mpi), INTENT(IN) :: n
2796
2797 REAL(mpd), INTENT(IN) :: cu(n)
2798 REAL(mpd), INTENT(IN) :: a(n,p)
2799 REAL(mpd), INTENT(IN) :: s((p*p+p)/2)
2800 REAL(mpd), INTENT(OUT) :: x(n+p)
2801 REAL(mpd), INTENT(IN) :: y(n+p)
2802
2803 REAL(mpd) :: dsum
2804
2805 DO i=1,n+p
2806 x(i)=y(i)
2807 END DO
2808 DO i=1,n
2809 x(i)=x(i)*cu(i) ! u =C^{-1/2} y
2810 DO j=1,p
2811 x(n+j)=x(n+j)-a(i,j)*x(i) ! d - K u
2812 END DO
2813 END DO
2814
2815 jj=0
2816 DO j=1,p ! Cholesky solution for v
2817 dsum=x(n+j)
2818 DO k=1,j-1
2819 dsum=dsum-s(k+jj)*x(n+k) ! H v = d - K u
2820 END DO
2821 x(n+j)=dsum ! -> v
2822 jj=jj+j
2823 END DO
2824
2825 DO j=p,1,-1 ! solution for lambda
2826 dsum=x(n+j)*s(jj)
2827 kk=jj
2828 DO k=j+1,p
2829 dsum=dsum+s(kk+j)*x(n+k) ! D H^T lambda = -v
2830 kk=kk+k
2831 END DO
2832 x(n+j)=-dsum ! -> lambda
2833 jj=jj-j
2834 END DO
2835
2836 DO i=1,n ! u - K^T lambda
2837 DO j=1,p
2838 x(i)=x(i)-a(i,j)*x(n+j)
2839 END DO
2840 END DO
2841 DO i=1,n
2842 x(i)=x(i)*cu(i) ! x = C^{-1/2} u
2843 END DO
2844
2845END SUBROUTINE presol
2846
2880
2881SUBROUTINE precons(p,n,b,nm,c,cu,a,l,s,nrkd)
2882 USE mpdef
2883
2884 IMPLICIT NONE
2885 INTEGER(mpi) :: i
2886 INTEGER(mpi) :: i0
2887 INTEGER(mpi) :: ib
2888 INTEGER(mpi) :: ii
2889 INTEGER(mpi) :: j
2890 INTEGER(mpi) :: jj
2891 INTEGER(mpi) :: jk
2892 INTEGER(mpi) :: k
2893 INTEGER(mpi) :: kk
2894 INTEGER(mpl) :: ij
2895 INTEGER(mpi) :: np
2896 INTEGER(mpi) :: p1
2897 INTEGER(mpi) :: p2
2898
2899 INTEGER(mpi), INTENT(IN) :: p
2900 INTEGER(mpi), INTENT(IN) :: n
2901 INTEGER(mpi), INTENT(IN) :: b
2902 INTEGER(mpl), INTENT(IN) :: nm
2903 REAL(mpd), INTENT(IN) :: c(n)
2904 REAL(mpd), INTENT(OUT) :: cu(n)
2905 REAL(mpd), INTENT(IN OUT) :: a(nm)
2906 INTEGER(mpi), INTENT(IN) :: l(p)
2907 REAL(mpd), INTENT(OUT) :: s((p*p+p)/2)
2908 INTEGER(mpi), INTENT(OUT) :: nrkd
2909
2910 REAL(mpd) :: div
2911 REAL(mpd) :: ratio
2912
2913 nrkd=0
2914 DO i=1,(p*p+p)/2
2915 s(i)=0.0_mpd
2916 END DO
2917 DO i=1,n
2918 div=c(i) ! copy
2919 IF (div > 0.0_mpd) THEN
2920 cu(i)=1.0_mpd/sqrt(div)
2921 ELSE
2922 cu(i)=0.0_mpd
2923 nrkd=nrkd+1
2924 END IF
2925 END DO
2926
2927 ij=0
2928 DO ib=1,b ! loop over blocks
2929 p1=l(ib*4-3) ! first constraint in block
2930 p2=l(ib*4-2) ! last constraint in block
2931 i0=l(ib*4-1) ! parameter offset in block
2932 np=l(ib*4) ! number of parameters in block
2933 DO i=1,np ! loop over parameters in block
2934 ij=ij+1
2935 DO j=p1,p2
2936 jk=(j*j-j)/2+p1
2937 a(ij+(j-p1)*np)=a(ij+(j-p1)*np)*cu(i+i0) ! K = A C^{-1/2}
2938 DO k=p1,j
2939 s(jk)=s(jk)+a(ij+(j-p1)*np)*a(ij+(k-p1)*np) ! S = symmetric matrix K K^T
2940 jk=jk+1
2941 END DO
2942 END DO
2943 END DO
2944 ij=ij+(p2-p1)*np
2945 END DO
2946
2947 ii=0
2948 DO i=1,p ! S -> H D H^T (Cholesky)
2949 ii=ii+i
2950 IF(s(ii) /= 0.0_mpd) s(ii)=1.0_mpd/s(ii)
2951 jj=ii
2952 DO j=i+1,p
2953 ratio=s(i+jj)*s(ii)
2954 kk=jj
2955 DO k=j,p
2956 s(kk+j)=s(kk+j)-s(kk+i)*ratio
2957 kk=kk+k
2958 END DO ! K
2959 s(i+jj)=ratio
2960 jj=jj+j
2961 END DO ! J
2962 END DO ! I
2963 RETURN
2964END SUBROUTINE precons
2965
2979
2980SUBROUTINE presols(p,n,b,nm,cu,a,l,s,x,y) ! solution
2981 USE mpdef
2982
2983 IMPLICIT NONE
2984 INTEGER(mpi) :: i
2985 INTEGER(mpi) :: i0
2986 INTEGER(mpi) :: ib
2987 INTEGER(mpi) :: j
2988 INTEGER(mpi) :: jj
2989 INTEGER(mpi) :: k
2990 INTEGER(mpi) :: kk
2991 INTEGER(mpl) :: ij
2992 INTEGER(mpi) :: np
2993 INTEGER(mpi) :: p1
2994 INTEGER(mpi) :: p2
2995
2996 INTEGER(mpi), INTENT(IN) :: p
2997 INTEGER(mpi), INTENT(IN) :: n
2998 INTEGER(mpi), INTENT(IN) :: b
2999 INTEGER(mpl), INTENT(IN) :: nm
3000
3001 REAL(mpd), INTENT(IN) :: cu(n)
3002 REAL(mpd), INTENT(IN) :: a(nm)
3003 INTEGER(mpi), INTENT(IN) :: l(p)
3004 REAL(mpd), INTENT(IN) :: s((p*p+p)/2)
3005 REAL(mpd), INTENT(OUT) :: x(n+p)
3006 REAL(mpd), INTENT(IN) :: y(n+p)
3007
3008 REAL(mpd) :: dsum
3009
3010 DO i=1,n+p
3011 x(i)=y(i)
3012 END DO
3013 DO i=1,n
3014 x(i)=x(i)*cu(i) ! u =C^{-1/2} y
3015 END DO
3016 ij=0
3017 DO ib=1,b ! loop over blocks
3018 p1=l(ib*4-3) ! first constraint in block
3019 p2=l(ib*4-2) ! last constraint in block
3020 i0=l(ib*4-1) ! parameter offset in block
3021 np=l(ib*4) ! number of parameters in block
3022 DO j=p1,p2 ! loop over constraints in block
3023 DO i=1,np ! loop over parameters in block
3024 ij=ij+1
3025 x(n+j)=x(n+j)-a(ij)*x(i+i0) ! d - K u
3026 END DO
3027 END DO
3028 END DO
3029
3030 jj=0
3031 DO j=1,p ! Cholesky solution for v
3032 dsum=x(n+j)
3033 DO k=1,j-1
3034 dsum=dsum-s(k+jj)*x(n+k) ! H v = d - K u
3035 END DO
3036 x(n+j)=dsum ! -> v
3037 jj=jj+j
3038 END DO
3039
3040 DO j=p,1,-1 ! solution for lambda
3041 dsum=x(n+j)*s(jj)
3042 kk=jj
3043 DO k=j+1,p
3044 dsum=dsum+s(kk+j)*x(n+k) ! D H^T lambda = -v
3045 kk=kk+k
3046 END DO
3047 x(n+j)=-dsum ! -> lambda
3048 jj=jj-j
3049 END DO
3050
3051 ij=0
3052 DO ib=1,b ! loop over blocks
3053 p1=l(ib*4-3) ! first constraint in block
3054 p2=l(ib*4-2) ! last constraint in block
3055 i0=l(ib*4-1) ! parameter offset in block
3056 np=l(ib*4) ! number of parameters in block
3057 DO j=p1,p2 ! loop over constraints in block
3058 DO i=1,np ! loop over parameters in block
3059 ij=ij+1
3060 x(i+i0)=x(i+i0)-a(ij)*x(n+j) ! u - K^T lambda
3061 END DO
3062 END DO
3063 END DO
3064
3065 DO i=1,n
3066 x(i)=x(i)*cu(i) ! x = C^{-1/2} u
3067 END DO
3068
3069END SUBROUTINE presols
3070
3071! 090817 C. Kleinwort, DESY-FH1
3116SUBROUTINE sqmibb(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
3117 USE mpdef
3118
3119 ! REAL(mpd) scratch arrays:
3120 ! VBND(N*(NBND+1)) = storage of band part
3121 ! VBDR(N* NBDR) = storage of border part
3122 ! AUX (N* NBDR) = intermediate results
3123
3124 ! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only)
3125
3126 IMPLICIT NONE
3127 INTEGER(mpi) :: i
3128 INTEGER(mpi) :: ib
3129 INTEGER(mpi) :: ij
3130 INTEGER(mpi) :: ioff
3131 INTEGER(mpi) :: ip
3132 INTEGER(mpi) :: ip1
3133 INTEGER(mpi) :: ip2
3134 INTEGER(mpi) :: is
3135 INTEGER(mpi) :: j
3136 INTEGER(mpi) :: j0
3137 INTEGER(mpi) :: jb
3138 INTEGER(mpi) :: joff
3139 INTEGER(mpi) :: mp1
3140 INTEGER(mpi) :: nb1
3141 INTEGER(mpi) :: nmb
3142 INTEGER(mpi) :: npri
3143 INTEGER(mpi) :: nrankb
3144
3145 INTEGER(mpi), INTENT(IN) :: n
3146 REAL(mpd), INTENT(IN OUT) :: v((n*n+n)/2)
3147 REAL(mpd), INTENT(OUT) :: b(n)
3148 INTEGER(mpi), INTENT(IN) :: nbdr
3149 INTEGER(mpi), INTENT(IN) :: nbnd
3150 INTEGER(mpi), INTENT(IN) :: inv
3151 INTEGER(mpi), INTENT(OUT) :: nrank
3152
3153 REAL(mpd), INTENT(OUT) :: vbnd(n*(nbnd+1))
3154 REAL(mpd), INTENT(OUT) :: vbdr(n*nbdr)
3155 REAL(mpd), INTENT(OUT) :: aux(n*nbdr)
3156 REAL(mpd), INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
3157 REAL(mpd), INTENT(OUT) :: vzru(nbdr)
3158 REAL(mpd), INTENT(OUT) :: scdiag(nbdr)
3159 INTEGER(mpi), INTENT(OUT) :: scflag(nbdr)
3160
3161 SAVE npri
3162 DATA npri / 100 /
3163 ! ...
3164 nrank=0
3165 nb1=nbdr+1
3166 mp1=nbnd+1
3167 nmb=n-nbdr
3168 ! copy band part
3169 DO i=nb1,n
3170 ip=(i*(i+1))/2
3171 is=0
3172 DO j=i,min(n,i+nbnd)
3173 ip=ip+is
3174 is=j
3175 ib=j-i+1
3176 vbnd(ib+(i-nb1)*mp1)=v(ip)
3177 END DO
3178 END DO
3179 ! copy border part
3180 IF (nbdr > 0) THEN
3181 ioff=0
3182 DO i=1,nbdr
3183 ip=(i*(i+1))/2
3184 is=0
3185 DO j=i,n
3186 ip=ip+is
3187 is=j
3188 vbdr(ioff+j)=v(ip)
3189 END DO
3190 ioff=ioff+n
3191 END DO
3192 END IF
3193
3194 CALL dbcdec(vbnd,mp1,nmb,aux)
3195 ! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable
3196 ! CALL DBCPRB(VBND,MP1,NMB)
3197 ip=1
3198 DO i=1, nmb
3199 IF (vbnd(ip) <= 0.0_mpd) THEN
3200 npri=npri-1
3201 IF (npri >= 0) THEN
3202 IF (vbnd(ip) == 0.0_mpd) THEN
3203 print *, ' SQMIBB matrix singular', n, nbdr, nbnd
3204 ELSE
3205 print *, ' SQMIBB matrix not positive definite', n, nbdr, nbnd
3206 END IF
3207 END IF
3208 ! return zeros
3209 DO ip=1,n
3210 b(ip)=0.0_mpd
3211 END DO
3212 DO ip=1,(n*n+n)/2
3213 v(ip)=0.0_mpd
3214 END DO
3215 RETURN
3216 END IF
3217 ip=ip+mp1
3218 END DO
3219 nrank=nmb
3220
3221 IF (nbdr == 0) THEN ! special case NBDR=0
3222
3223 CALL dbcslv(vbnd,mp1,nmb,b,b)
3224 IF (inv > 0) THEN
3225 IF (inv > 1) THEN
3226 CALL dbcinv(vbnd,mp1,nmb,v)
3227 ELSE
3228 CALL dbcinb(vbnd,mp1,nmb,v)
3229 END IF
3230 END IF
3231
3232 ELSE ! general case NBDR>0
3233
3234 ioff=nb1
3235 DO ib=1,nbdr
3236 ! solve for aux. vectors
3237 CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff),aux(ioff))
3238 ! zT ru
3239 vzru(ib)=b(ib)
3240 DO i=0,nmb-1
3241 vzru(ib)=vzru(ib)-b(nb1+i)*aux(ioff+i)
3242 END DO
3243 ioff=ioff+n
3244 END DO
3245 ! solve for band part only
3246 CALL dbcslv(vbnd,mp1,nmb,b(nb1),b(nb1))
3247 ! Ck - cT z
3248 ip=0
3249 ioff=nb1
3250 DO ib=1,nbdr
3251 joff=nb1
3252 DO jb=1,ib
3253 ip=ip+1
3254 vbk(ip)=v(ip)
3255 DO i=0,nmb-1
3256 vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
3257 END DO
3258 joff=joff+n
3259 END DO
3260 ioff=ioff+n
3261 END DO
3262 ! solve border part
3263 CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
3264 IF (nrankb == nbdr) THEN
3265 nrank=nrank+nbdr
3266 ELSE
3267 npri=npri-1
3268 IF (npri >= 0) print *, ' SQMIBB undef border ', n, nbdr, nbnd, nrankb
3269 DO ib=1,nbdr
3270 vzru(ib)=0.0_mpd
3271 END DO
3272 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
3273 vbk(ip)=0.0_mpd
3274 END DO
3275 END IF
3276 ! smoothed data points
3277 ioff=nb1
3278 DO ib=1, nbdr
3279 b(ib) = vzru(ib)
3280 DO i=0,nmb-1
3281 b(nb1+i)=b(nb1+i)-b(ib)*aux(ioff+i)
3282 END DO
3283 ioff=ioff+n
3284 END DO
3285 ! inverse requested ?
3286 IF (inv > 0) THEN
3287 IF (inv > 1) THEN
3288 CALL dbcinv(vbnd,mp1,nmb,v)
3289 ELSE
3290 CALL dbcinb(vbnd,mp1,nmb,v)
3291 END IF
3292 ! expand/correct from NMB to N
3293 ip1=(nmb*nmb+nmb)/2
3294 ip2=(n*n+n)/2
3295 DO i=nmb-1,0,-1
3296 j0=0
3297 IF (inv == 1) j0=max(0,i-nbnd)
3298 DO j=i,j0,-1
3299 v(ip2)=v(ip1)
3300 ioff=nb1
3301 DO ib=1,nbdr
3302 joff=nb1
3303 DO jb=1,nbdr
3304 ij=max(ib,jb)
3305 ij=(ij*ij-ij)/2+min(ib,jb)
3306 v(ip2)=v(ip2)+vbk(ij)*aux(ioff+i)*aux(joff+j)
3307 joff=joff+n
3308 END DO
3309 ioff=ioff+n
3310 END DO
3311 ip1=ip1-1
3312 ip2=ip2-1
3313 END DO
3314 ip1=ip1-j0
3315 ip2=ip2-j0
3316
3317 DO ib=nbdr,1,-1
3318 v(ip2)=0.0_mpd
3319 joff=nb1
3320 DO jb=1,nbdr
3321 ij=max(ib,jb)
3322 ij=(ij*ij-ij)/2+min(ib,jb)
3323 v(ip2)=v(ip2)-vbk(ij)*aux(i+joff)
3324 joff=joff+n
3325 END DO
3326 ip2=ip2-1
3327 END DO
3328 END DO
3329
3330 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
3331 v(ip2)=vbk(ip)
3332 ip2=ip2-1
3333 END DO
3334
3335 END IF
3336 END IF
3337
3338END SUBROUTINE sqmibb
3339
3340! 181105 C. Kleinwort, DESY-BELLE
3372SUBROUTINE sqmibb2(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
3373 USE mpdef
3374
3375 ! REAL(mpd) scratch arrays:
3376 ! VBND(N*(NBND+1)) = storage of band part
3377 ! VBDR(N* NBDR) = storage of border part
3378 ! AUX (N* NBDR) = intermediate results
3379
3380 ! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only)
3381
3382 IMPLICIT NONE
3383 INTEGER(mpi) :: i
3384 INTEGER(mpi) :: ib
3385 INTEGER(mpi) :: ij
3386 INTEGER(mpi) :: ioff
3387 INTEGER(mpi) :: ip
3388 INTEGER(mpi) :: ip1
3389 INTEGER(mpi) :: is
3390 INTEGER(mpi) :: j
3391 INTEGER(mpi) :: j0
3392 INTEGER(mpi) :: jb
3393 INTEGER(mpi) :: joff
3394 INTEGER(mpi) :: koff
3395 INTEGER(mpi) :: mp1
3396 INTEGER(mpi) :: nb1
3397 INTEGER(mpi) :: nmb
3398 INTEGER(mpi) :: npri
3399 INTEGER(mpi) :: nrankb
3400
3401 INTEGER(mpi), INTENT(IN) :: n
3402 REAL(mpd), INTENT(IN OUT) :: v((n*n+n)/2)
3403 REAL(mpd), INTENT(OUT) :: b(n)
3404 INTEGER(mpi), INTENT(IN) :: nbdr
3405 INTEGER(mpi), INTENT(IN) :: nbnd
3406 INTEGER(mpi), INTENT(IN) :: inv
3407 INTEGER(mpi), INTENT(OUT) :: nrank
3408
3409 REAL(mpd), INTENT(OUT) :: vbnd(n*(nbnd+1))
3410 REAL(mpd), INTENT(OUT) :: vbdr(n*nbdr)
3411 REAL(mpd), INTENT(OUT) :: aux(n*nbdr)
3412 REAL(mpd), INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
3413 REAL(mpd), INTENT(OUT) :: vzru(nbdr)
3414 REAL(mpd), INTENT(OUT) :: scdiag(nbdr)
3415 INTEGER(mpi), INTENT(OUT) :: scflag(nbdr)
3416
3417 SAVE npri
3418 DATA npri / 100 /
3419 ! ...
3420 nrank=0
3421 mp1=nbnd+1
3422 nmb=n-nbdr
3423 nb1=nmb+1
3424 ! copy band part
3425 DO i=1,nmb
3426 ip=(i*(i+1))/2
3427 is=0
3428 DO j=i,min(nmb,i+nbnd)
3429 ip=ip+is
3430 is=j
3431 ib=j-i+1
3432 vbnd(ib+(i-1)*mp1)=v(ip)
3433 END DO
3434 END DO
3435 ! copy border part
3436 IF (nbdr > 0) THEN
3437 ioff=0
3438 DO i=nb1,n
3439 ip=(i*(i-1))/2
3440 DO j=1,i
3441 vbdr(ioff+j)=v(ip+j)
3442 END DO
3443 ioff=ioff+n
3444 END DO
3445 END IF
3446
3447 CALL dbcdec(vbnd,mp1,nmb,aux)
3448 ! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable
3449 ! CALL DBCPRB(VBND,MP1,NMB)
3450 ip=1
3451 DO i=1, nmb
3452 IF (vbnd(ip) <= 0.0_mpd) THEN
3453 npri=npri-1
3454 IF (npri >= 0) THEN
3455 IF (vbnd(ip) == 0.0_mpd) THEN
3456 print *, ' SQMIBB2 matrix singular', n, nbdr, nbnd
3457 ELSE
3458 print *, ' SQMIBB2 matrix not positive definite', n, nbdr, nbnd
3459 END IF
3460 END IF
3461 ! return zeros
3462 DO ip=1,n
3463 b(ip)=0.0_mpd
3464 END DO
3465 DO ip=1,(n*n+n)/2
3466 v(ip)=0.0_mpd
3467 END DO
3468 RETURN
3469 END IF
3470 ip=ip+mp1
3471 END DO
3472 nrank=nmb
3473
3474 IF (nbdr == 0) THEN ! special case NBDR=0
3475
3476 CALL dbcslv(vbnd,mp1,nmb,b,b)
3477 IF (inv > 0) THEN
3478 IF (inv > 1) THEN
3479 CALL dbcinv(vbnd,mp1,nmb,v)
3480 ELSE
3481 CALL dbcinb(vbnd,mp1,nmb,v)
3482 END IF
3483 END IF
3484
3485 ELSE ! general case NBDR>0
3486
3487 ioff=0
3488 DO ib=1,nbdr
3489 ! solve for aux. vectors
3490 CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff+1),aux(ioff+1))
3491 ! zT ru
3492 vzru(ib)=b(nmb+ib)
3493 DO i=1,nmb
3494 vzru(ib)=vzru(ib)-b(i)*aux(ioff+i)
3495 END DO
3496 ioff=ioff+n
3497 END DO
3498 ! solve for band part only
3499 CALL dbcslv(vbnd,mp1,nmb,b,b)
3500 ! Ck - cT z
3501 ip=0
3502 ioff=0
3503 koff=nmb
3504 DO ib=1,nbdr
3505 joff=0
3506 DO jb=1,ib
3507 ip=ip+1
3508 vbk(ip)=vbdr(koff+jb)
3509 DO i=1,nmb
3510 vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
3511 END DO
3512 joff=joff+n
3513 END DO
3514 ioff=ioff+n
3515 koff=koff+n
3516 END DO
3517
3518 ! solve border part
3519 CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
3520 IF (nrankb == nbdr) THEN
3521 nrank=nrank+nbdr
3522 ELSE
3523 npri=npri-1
3524 IF (npri >= 0) print *, ' SQMIBB2 undef border ', n, nbdr, nbnd, nrankb
3525 DO ib=1,nbdr
3526 vzru(ib)=0.0_mpd
3527 END DO
3528 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
3529 vbk(ip)=0.0_mpd
3530 END DO
3531 END IF
3532 ! smoothed data points
3533 ioff=0
3534 DO ib=1, nbdr
3535 DO i=1,nmb
3536 b(i)=b(i)-vzru(ib)*aux(ioff+i)
3537 END DO
3538 ioff=ioff+n
3539 b(nmb+ib)=vzru(ib)
3540 END DO
3541 ! inverse requested ?
3542 IF (inv > 0) THEN
3543 IF (inv > 1) THEN
3544 CALL dbcinv(vbnd,mp1,nmb,v)
3545 ELSE
3546 CALL dbcinb(vbnd,mp1,nmb,v)
3547 END IF
3548 ! assemble band and border
3549 IF (nbdr > 0) THEN
3550 ! band part
3551 ip1=(nmb*nmb+nmb)/2
3552 DO i=nmb-1,0,-1
3553 j0=0
3554 IF (inv == 1) j0=max(0,i-nbnd)
3555 DO j=i,j0,-1
3556 ioff=1
3557 DO ib=1,nbdr
3558 joff=1
3559 DO jb=1,nbdr
3560 ij=max(ib,jb)
3561 ij=(ij*ij-ij)/2+min(ib,jb)
3562 v(ip1)=v(ip1)+vbk(ij)*aux(ioff+i)*aux(joff+j)
3563 joff=joff+n
3564 END DO
3565 ioff=ioff+n
3566 END DO
3567 ip1=ip1-1
3568 END DO
3569 ip1=ip1-j0
3570 END DO
3571 ! border part
3572 ip1=(nmb*nmb+nmb)/2
3573 ip=0
3574 DO ib=1,nbdr
3575 DO i=1,nmb
3576 ip1=ip1+1
3577 v(ip1)=0.0_mpd
3578 joff=0
3579 DO jb=1,nbdr
3580 ij=max(ib,jb)
3581 ij=(ij*ij-ij)/2+min(ib,jb)
3582 v(ip1)=v(ip1)-vbk(ij)*aux(i+joff)
3583 joff=joff+n
3584 END DO
3585 END DO
3586 DO jb=1,ib
3587 ip1=ip1+1
3588 ip=ip+1
3589 v(ip1)=vbk(ip)
3590 END DO
3591 END DO
3592
3593 END IF
3594 END IF
3595 END IF
3596
3597END SUBROUTINE sqmibb2
dbcslv
subroutine dbcslv(w, mp1, n, b, x)
solution B -> X
Definition: Dbandmatrix.f90:504
dbcinv
subroutine dbcinv(w, mp1, n, vs)
VS = inverse symmetric matrix.
Definition: Dbandmatrix.f90:614
dbcinb
subroutine dbcinb(w, mp1, n, vs)
VS = band part of inverse symmetric matrix.
Definition: Dbandmatrix.f90:577
dbcdec
subroutine dbcdec(w, mp1, n, aux)
Decomposition of symmetric band matrix.
Definition: Dbandmatrix.f90:460
monpgs
subroutine monpgs(i)
Progress monitoring.
Definition: mpmon.f90:68
cholin
subroutine cholin(g, v, n)
Inversion after decomposition.
Definition: mpnum.f90:838
sqmibb2
subroutine sqmibb2(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Band bordered matrix.
Definition: mpnum.f90:3373
vabmmm
subroutine vabmmm(n, val, ilptr)
Variable band matrix print minimum and maximum.
Definition: mpnum.f90:1125
dbavat
subroutine dbavat(v, a, w, n, m, iopt)
A V AT product (similarity).
Definition: mpnum.f90:1390
sqminl
subroutine sqminl(v, b, n, nrank, diag, next, vk, mon)
Matrix inversion for LARGE matrices.
Definition: mpnum.f90:231
precon
subroutine precon(p, n, c, cu, a, s, nrkd)
Constrained preconditioner, decomposition.
Definition: mpnum.f90:2711
devsol
subroutine devsol(n, diag, u, b, x, work)
Solution by diagonalization.
Definition: mpnum.f90:650
dbsvxl
subroutine dbsvxl(v, a, b, n)
Product LARGE symmetric matrix, vector.
Definition: mpnum.f90:1309
devrot
subroutine devrot(n, diag, u, v, work, iwork)
Diagonalization.
Definition: mpnum.f90:370
sort22l
subroutine sort22l(a, b, n)
Quick sort 2 with index.
Definition: mpnum.f90:1982
dbmprv
subroutine dbmprv(lun, x, v, n)
Print symmetric matrix, vector.
Definition: mpnum.f90:1541
dbaxpy
subroutine dbaxpy(n, da, dx, dy)
Multiply, addition.
Definition: mpnum.f90:1234
dbavats
subroutine dbavats(v, a, is, w, n, m, iopt, sc)
A V AT product (similarity, sparse).
Definition: mpnum.f90:1471
equslv
subroutine equslv(n, m, c, india, x)
Solution of equilibrium systems (after decomposition).
Definition: mpnum.f90:2423
heapf
subroutine heapf(a, n)
Heap sort direct (real).
Definition: mpnum.f90:1660
dbdot
real(mpd) function dbdot(n, dx, dy)
Dot product.
Definition: mpnum.f90:1203
equdec
subroutine equdec(n, m, ls, c, india, nrkd, nrkd2)
Decomposition of equilibrium systems.
Definition: mpnum.f90:2356
vabdec
subroutine vabdec(n, val, ilptr)
Variable band matrix decomposition.
Definition: mpnum.f90:1013
chslv2
subroutine chslv2(g, x, n)
Solve A*x=b using Cholesky decomposition.
Definition: mpnum.f90:954
vabslv
subroutine vabslv(n, val, ilptr, x)
Variable band matrix solution.
Definition: mpnum.f90:1162
choldc
subroutine choldc(g, n)
Cholesky decomposition.
Definition: mpnum.f90:746
sort1k
subroutine sort1k(a, n)
Quick sort 1.
Definition: mpnum.f90:1715
sqminv
subroutine sqminv(v, b, n, nrank, diag, next)
Matrix inversion and solution.
Definition: mpnum.f90:98
lltbwd
subroutine lltbwd(n, c, india, x)
Backward solution.
Definition: mpnum.f90:2316
presols
subroutine presols(p, n, b, nm, cu, a, l, s, x, y)
Constrained (sparse) preconditioner, solution.
Definition: mpnum.f90:2981
dbgax
subroutine dbgax(a, x, y, m, n)
Multiply general M-by-N matrix A and N-vector X.
Definition: mpnum.f90:1352
sqmibb
subroutine sqmibb(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Bordered band matrix.
Definition: mpnum.f90:3117
lltfwds
subroutine lltfwds(n, c, india, x, i0, ns)
Forward solution (sparse).
Definition: mpnum.f90:2276
presol
subroutine presol(p, n, cu, a, s, x, y)
Constrained preconditioner, solution.
Definition: mpnum.f90:2785
dbprv
subroutine dbprv(lun, v, n)
Print symmetric matrix.
Definition: mpnum.f90:1616
devinv
subroutine devinv(n, diag, u, v)
Inversion by diagonalization.
Definition: mpnum.f90:697
lltdec
subroutine lltdec(n, c, india, nrkd, iopt)
LLT decomposition.
Definition: mpnum.f90:2156
equdecs
subroutine equdecs(n, m, b, nm, ls, c, india, l, nrkd, nrkd2)
Decomposition of (sparse) equilibrium systems.
Definition: mpnum.f90:2487
cholsl
subroutine cholsl(g, x, n)
Solution after decomposition.
Definition: mpnum.f90:792
chdec2
subroutine chdec2(g, n, nrank, evmax, evmin, mon)
Cholesky decomposition (LARGE pos.
Definition: mpnum.f90:892
sort2k
subroutine sort2k(a, n)
Quick sort 2.
Definition: mpnum.f90:1800
devsig
subroutine devsig(n, diag, u, b, coef)
Calculate significances.
Definition: mpnum.f90:612
chindl
real(mps) function chindl(n, nd)
Chi2/ndf cuts.
Definition: mpnum.f90:2079
lltfwd
subroutine lltfwd(n, c, india, x)
Forward solution.
Definition: mpnum.f90:2238
dbsvx
subroutine dbsvx(v, a, b, n)
Product symmetric matrix, vector.
Definition: mpnum.f90:1265
equslvs
subroutine equslvs(n, m, b, nm, c, india, l, x)
Solution of (sparse) equilibrium systems (after decomposition).
Definition: mpnum.f90:2614
precons
subroutine precons(p, n, b, nm, c, cu, a, l, s, nrkd)
Constrained (sparse) preconditioner, decomposition.
Definition: mpnum.f90:2882
sort2i
subroutine sort2i(a, n)
Quick sort 2 with index.
Definition: mpnum.f90:1893
mpdef
Definition of constants.
Definition: mpdef.f90:24
mpdef::mpd
integer, parameter mpd
double precision
Definition: mpdef.f90:38
mpdef::mpl
integer, parameter mpl
long integer
Definition: mpdef.f90:36
mpdef::mps
integer, parameter mps
single precision
Definition: mpdef.f90:37
mpdef::mpi
integer, parameter mpi
integer
Definition: mpdef.f90:35
peend
subroutine peend(icode, cmessage)
Print exit code.
Definition: pede.f90:12984
Generated on Fri Apr 12 2024 12:19:37 for Millepede-II by doxygen 1.9.3