25integer :: m_pds_neqns1
26integer :: m_pds_neqns2
27integer :: m_pds_neqnsd
29real(8),
pointer :: m_pds_d(:,:)
30type(ccls_matrix),
pointer :: m_pds_c1, m_pds_c2
31integer,
pointer :: m_pds_part(:)
32logical :: m_pds_isready_dc = .false.
35integer,
parameter :: m_pds_lenv=80000000
36real(8),
pointer :: m_pds_v(:)
37logical :: m_pds_isready_v = .false.
59character(132),
parameter::
version=
'sp_LINEQ: $Revision: 1.2 $ $Date: 2007/11/21 02:13:45 $'
62type (sp_matrix) :: sp_mat
74 call sp_direct_root(sp_mat)
77 call sp_direct_trunk()
91subroutine sp_direct_root(sp_mat)
99character(132),
parameter ::
version=
'sp_direct_root: $Revision: 1.2 $ $Date: 2007/11/21 02:13:45 $'
104type (sp_matrix) :: sp_mat
111real(8),
allocatable,
dimension(:) :: r
114integer :: nprof, nprof2
119integer,
allocatable,
dimension(:) :: iperm
126real(8),
allocatable,
dimension(:,:) :: d1, d2
131integer :: i,j,k,l,m,n
134integer,
dimension(MPI_STATUS_SIZE) :: istatus
138logical,
parameter :: elap=.true.
141real(8) :: t20s, t20e, t21s, t21e
142real(8) :: t30s, t30e, t31s, t31e
148integer,
allocatable,
dimension(:),
target :: part
149real(8),
allocatable,
dimension(:,:),
target :: d
150integer :: neqns1, neqns2, neqnsd, nd
151real(8),
allocatable,
dimension(:) :: oldb
159call elapout(
'sp_direct_root: begin')
170nprof=sp_mat%ipoi(sp_mat%ncol+1)-1
171nprof2=nprof/2+a0%neqns
172allocate(a0%irow(nprof/2+a0%neqns),a0%jcol(nprof/2+a0%neqns),a0%val(nndeg,nprof/2+a0%neqns))
178 a0%val(1:nndeg,ll)=sp_mat%d(jj,1:nndeg)
181 ii=sp_mat%irowno(loc)
187 a0%val(1:nndeg,ll)=sp_mat%elm(loc,1:nndeg)
195allocate(r(a0%neqns*ndeg))
198 r(ndeg*(i-1)+j)=sp_mat%b(j,i)
203call elapout(
'sp_direct_root: get matrix information done ')
208call elapout(
'sp_direct_root: begin divide matrix ')
211allocate(part(a0%neqns))
212allocate(iperm(a0%neqns))
215call mkpart(a0%neqns, a0%nttbr, ndeg, a0%irow, a0%jcol, neqns1, neqns2, neqnsd, part, iperm)
221call divmat(a0, part, iperm, neqns1, neqns2, neqnsd, ndeg, a1, a2, c1, c2, d)
223call elapout(
'sp_direct_root: end divide matrix')
233call elapout(
'sp_direct_root: begin send a1 ')
235call mpi_send(a1%neqns, 1, mpi_integer,ip1,1,mpi_comm_world,ierr)
236call mpi_send(a1%nttbr, 1, mpi_integer,ip1,1,mpi_comm_world,ierr)
237call mpi_send(a1%ndeg, 1, mpi_integer,ip1,1,mpi_comm_world,ierr)
239call mpi_send(a1%irow, a1%nttbr, mpi_integer, ip1,1,mpi_comm_world,ierr)
240call mpi_send(a1%jcol, a1%nttbr, mpi_integer, ip1,1,mpi_comm_world,ierr)
241call mpi_send(a1%val, a1%nttbr*nndeg, mpi_real8, ip1,1,mpi_comm_world,ierr)
242call elapout(
'sp_direct_root: end send a1 ')
244call elapout(
'sp_direct_root: begin send a2 ')
246call mpi_send(a2%neqns, 1, mpi_integer,ip2,1,mpi_comm_world,ierr)
247call mpi_send(a2%nttbr, 1, mpi_integer,ip2,1,mpi_comm_world,ierr)
248call mpi_send(a2%ndeg, 1, mpi_integer,ip2,1,mpi_comm_world,ierr)
250call mpi_send(a2%irow, a2%nttbr, mpi_integer, ip2,1,mpi_comm_world,ierr)
251call mpi_send(a2%jcol, a2%nttbr, mpi_integer, ip2,1,mpi_comm_world,ierr)
252call mpi_send(a2%val, a2%nttbr*nndeg, mpi_real8, ip2,1,mpi_comm_world,ierr)
253call elapout(
'sp_direct_root: end send a2 ')
266call elapout(
'sp_direct_root: begin send c1 ')
267call mpi_send(c1%neqns, 1, mpi_integer, ip1,1,mpi_comm_world,ierr)
268call mpi_send(c1%ncol, 1, mpi_integer, ip1,1,mpi_comm_world,ierr)
269call mpi_send(c1%nttbr, 1, mpi_integer, ip1,1,mpi_comm_world,ierr)
270call mpi_send(c1%ndeg, 1, mpi_integer, ip1,1,mpi_comm_world,ierr)
271call mpi_send(c1%ia, c1%ncol+1, mpi_integer, ip1,1,mpi_comm_world,ierr)
272call mpi_send(c1%ja, c1%nttbr, mpi_integer, ip1,1,mpi_comm_world,ierr)
273call mpi_send(c1%val, c1%nttbr*nndeg, mpi_real8, ip1,1,mpi_comm_world,ierr)
274call elapout(
'sp_direct_root: end send c1 ')
276call elapout(
'sp_direct_root: begin send c2')
277call mpi_send(c2%neqns, 1, mpi_integer, ip2,1,mpi_comm_world,ierr)
278call mpi_send(c2%ncol, 1, mpi_integer, ip2,1,mpi_comm_world,ierr)
279call mpi_send(c2%nttbr, 1, mpi_integer, ip2,1,mpi_comm_world,ierr)
280call mpi_send(c2%ndeg, 1, mpi_integer, ip2,1,mpi_comm_world,ierr)
281call mpi_send(c2%ia, c2%ncol+1, mpi_integer, ip2,1,mpi_comm_world,ierr)
282call mpi_send(c2%ja, c2%nttbr, mpi_integer, ip2,1,mpi_comm_world,ierr)
283call mpi_send(c2%val, c2%nttbr*nndeg, mpi_real8, ip2,1,mpi_comm_world,ierr)
284call elapout(
'sp_direct_root: end send c2')
286allocate(d1(nd,nd), d2(nd,nd))
290call elapout(
'sp_direct_root: wait until receive D1, D2')
291call mpi_recv(d1, nd*nd, mpi_real8, ip1, 1,mpi_comm_world, istatus, ierr)
292call mpi_recv(d2, nd*nd, mpi_real8, ip2, 1,mpi_comm_world, istatus, ierr)
293call elapout(
'sp_direct_root: receive d1, d2 from children')
296call elapout(
'sp_direct_root: modify D and LDU decompose it.')
303call elapout(
'sp_direct_root: LDU decompose for denssol done.')
306m_pds_neqns = a0%neqns
317m_pds_isready_dc=.true.
329allocate(oldb(a0%neqns*ndeg))
331call elapout(
'sp_direct_root: entering trunksolve')
333call elapout(
'sp_direct_root: end trunksolve')
335call verif0(a0%neqns, a0%ndeg, a0%nttbr, a0%irow, a0%jcol, a0%val, oldb, r)
340 sp_mat%x(j,i)=r(ndeg*(i-1)+j)
346call elapout(
'sp_direct_root: send stop flag to children')
348call mpi_send(do_calc, 1, mpi_logical, ip1,1,mpi_comm_world,ierr)
349call mpi_send(do_calc, 1, mpi_logical, ip2,1,mpi_comm_world,ierr)
351call elapout(
'sp_direct_root: end')
363end subroutine sp_direct_root
367subroutine sp_direct_trunk()
375character(132),
parameter ::
version=
'sp_direct_trunk: $Revision: 1.2 $ $Date: 2007/11/21 02:13:45 $'
387integer,
allocatable,
dimension(:) :: iperm
394real(8),
allocatable,
dimension(:,:) :: d1, d2
399integer :: i,j,k,l,m,n
402integer,
dimension(MPI_STATUS_SIZE) :: istatus
405logical,
parameter :: elap=.true.
408real(8) :: t20s, t20e, t21s, t21e
409real(8) :: t30s, t30e, t31s, t31e
415integer,
allocatable,
dimension(:),
target :: part
416real(8),
allocatable,
dimension(:,:),
target :: d
417integer :: neqns1, neqns2, neqnsd, nd
424call elapout(
'sp_direct_trunk: begin')
425call elapout(
'sp_direct_trunk: waiting matrix via MPI')
427call mpi_recv(a0%neqns, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
428call mpi_recv(a0%nttbr, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
429call mpi_recv(a0%ndeg, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
430allocate(a0%irow(a0%nttbr))
431allocate(a0%jcol(a0%nttbr))
432allocate(a0%val(a0%ndeg*a0%ndeg,a0%nttbr))
434call elapout(
'sp_direct_trunk: begen get matrix via MPI')
435call mpi_recv(a0%irow, a0%nttbr, mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
436call mpi_recv(a0%jcol, a0%nttbr, mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
437call mpi_recv(a0%val, a0%nttbr*a0%ndeg*a0%ndeg, mpi_real8,imp,1,mpi_comm_world,istatus,ierr)
439call elapout(
'sp_direct_trunk: end get matrix via MPI')
448call elapout(
'sp_direct_trunk: begin divide matrix')
450allocate(part(a0%neqns))
451allocate(iperm(a0%neqns))
454call mkpart(a0%neqns, a0%nttbr, ndeg, a0%irow, a0%jcol, neqns1, neqns2, neqnsd, part, iperm)
460call divmat(a0, part, iperm, neqns1, neqns2, neqnsd, ndeg, a1, a2, c1, c2, d)
461call elapout(
'sp_direct_trunk: end divide matrix')
469call elapout(
'sp_direct_trunk: begin send a1 ')
471call mpi_send(a1%neqns, 1, mpi_integer,ip1,1,mpi_comm_world,ierr)
472call mpi_send(a1%nttbr, 1, mpi_integer,ip1,1,mpi_comm_world,ierr)
473call mpi_send(a1%ndeg, 1, mpi_integer,ip1,1,mpi_comm_world,ierr)
475call mpi_send(a1%irow, a1%nttbr, mpi_integer, ip1,1,mpi_comm_world,ierr)
476call mpi_send(a1%jcol, a1%nttbr, mpi_integer, ip1,1,mpi_comm_world,ierr)
477call mpi_send(a1%val, a1%nttbr*nndeg, mpi_real8, ip1,1,mpi_comm_world,ierr)
478call elapout(
'sp_direct_trunk: end send a1 ')
480call elapout(
'sp_direct_trunk: begin send a2 ')
482call mpi_send(a2%neqns, 1, mpi_integer,ip2,1,mpi_comm_world,ierr)
483call mpi_send(a2%nttbr, 1, mpi_integer,ip2,1,mpi_comm_world,ierr)
484call mpi_send(a2%ndeg, 1, mpi_integer,ip2,1,mpi_comm_world,ierr)
486call mpi_send(a2%irow, a2%nttbr, mpi_integer, ip2,1,mpi_comm_world,ierr)
487call mpi_send(a2%jcol, a2%nttbr, mpi_integer, ip2,1,mpi_comm_world,ierr)
488call mpi_send(a2%val, a2%nttbr*nndeg, mpi_real8, ip2,1,mpi_comm_world,ierr)
489call elapout(
'sp_direct_trunk: end send a2')
500call elapout(
'sp_direct_trunk: send c1')
501call mpi_send(c1%neqns, 1, mpi_integer, ip1,1,mpi_comm_world,ierr)
502call mpi_send(c1%ncol, 1, mpi_integer, ip1,1,mpi_comm_world,ierr)
503call mpi_send(c1%nttbr, 1, mpi_integer, ip1,1,mpi_comm_world,ierr)
504call mpi_send(c1%ndeg, 1, mpi_integer, ip1,1,mpi_comm_world,ierr)
505call mpi_send(c1%ia, c1%ncol+1, mpi_integer, ip1,1,mpi_comm_world,ierr)
506call mpi_send(c1%ja, c1%nttbr, mpi_integer, ip1,1,mpi_comm_world,ierr)
507call mpi_send(c1%val, c1%nttbr*nndeg, mpi_real8, ip1,1,mpi_comm_world,ierr)
509call elapout(
'sp_direct_trunk: send c2')
510call mpi_send(c2%neqns, 1, mpi_integer, ip2,1,mpi_comm_world,ierr)
511call mpi_send(c2%ncol, 1, mpi_integer, ip2,1,mpi_comm_world,ierr)
512call mpi_send(c2%nttbr, 1, mpi_integer, ip2,1,mpi_comm_world,ierr)
513call mpi_send(c2%ndeg, 1, mpi_integer, ip2,1,mpi_comm_world,ierr)
514call mpi_send(c2%ia, c2%ncol+1, mpi_integer, ip2,1,mpi_comm_world,ierr)
515call mpi_send(c2%ja, c2%nttbr, mpi_integer, ip2,1,mpi_comm_world,ierr)
516call mpi_send(c2%val, c2%nttbr*nndeg, mpi_real8, ip2,1,mpi_comm_world,ierr)
518allocate(d1(nd,nd), d2(nd,nd))
522call elapout(
'sp_direct_trunk: wait until receive D1, D2')
523call mpi_recv(d1, nd*nd, mpi_real8, ip1, 1,mpi_comm_world, istatus, ierr)
524call mpi_recv(d2, nd*nd, mpi_real8, ip2, 1,mpi_comm_world, istatus, ierr)
525call elapout(
'sp_direct_trunk: receive d1, d2 from children')
528call elapout(
'sp_direct_trunk: modify D and LDU decompose it.')
535call elapout(
'sp_direct_trunk: LDU decompose for denssol done.')
538m_pds_neqns = a0%neqns
549m_pds_isready_dc=.true.
556call elapout(
'sp_direct_trunk: begin calcCtAC()')
558call elapout(
'sp_direct_trunk: end calcCtAC()')
565call elapout(
'sp_direct_trunk: begin recsolve()')
567call elapout(
'sp_direct_trunk: end recsolve()')
571end subroutine sp_direct_trunk
575subroutine sp_direct_leaf()
589integer,
dimension(MPI_STATUS_SIZE) :: istatus
598real(8) :: t20s, t20e, t21s, t21e
599real(8) :: t30s, t30e, t31s, t31e
609call elapout(
'sp_direct_leaf: begin')
610call elapout(
'sp_direct_leaf: waiting matrix via MPI')
612call mpi_recv(a%neqns, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
613call mpi_recv(a%nttbr, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
614call mpi_recv(a%ndeg, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
615allocate(a%irow(a%nttbr))
616allocate(a%jcol(a%nttbr))
617allocate(a%val(a%ndeg*a%ndeg,a%nttbr))
619call elapout(
'sp_direct_leaf: begin receive matrix')
621call mpi_recv(a%irow, a%nttbr, mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
622call mpi_recv(a%jcol, a%nttbr, mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
623call mpi_recv(a%val, a%nttbr*a%ndeg*a%ndeg, mpi_real8,imp,1,mpi_comm_world,istatus,ierr)
625call elapout(
'sp_direct_leaf: end get matrix via MPI')
633allocate(m_pds_v(m_pds_lenv))
635call elapout(
'sp_direct_leaf: begin matini()')
636call matini(a%neqns, a%nttbr, a%irow, a%jcol, m_pds_lenv, m_pds_v, ierr)
637if (ierr .ne. 0)
call errtrp(
'sp_direct_leaf: matini')
638call elapout(
'sp_direct_leaf: end matini()')
640call elapout(
'sp_direct_leaf: begin staij1()')
642 call staij1(0, a%irow(i), a%jcol(i), a%val(1,i), m_pds_v, a%ndeg, ierr)
643 if (ierr .ne. 0)
call errtrp(
'sp_direct_leaf: staij1')
645call elapout(
'sp_direct_leaf: end staij1()')
647call elapout(
'sp_direct_leaf: begin nufct0()')
648call nufct0(m_pds_v, ierr)
649if (ierr .ne. 0)
call errtrp(
'sp_direct_leaf: nufct0')
650call elapout(
'sp_direct_leaf: end nufct0()')
655m_pds_isready_v=.true.
663call elapout(
'sp_direct_leaf: entering calcCtAC')
667call elapout(
'sp_direct_leaf: exit calcCtAC')
674call elapout(
'sp_direct_leaf: entering recsolve')
678call elapout(
'sp_direct_leaf: exit recsolve')
691end subroutine sp_direct_leaf
705real(8),
dimension(m_pds_ndeg*m_pds_neqns) :: b
709integer,
dimension(MPI_STATUS_SIZE) :: istatus
714 call mpi_recv(do_calc, 1, mpi_logical, imp, 1,mpi_comm_world, istatus, ierr)
715 if (.false.==do_calc)
then
716 if (.true. == istrunk)
then
718 call mpi_send(do_calc, 1, mpi_logical, ip1,1,mpi_comm_world,ierr)
719 call mpi_send(do_calc, 1, mpi_logical, ip2,1,mpi_comm_world,ierr)
729 call mpi_recv(b, m_pds_ndeg*m_pds_neqns, mpi_real8, imp, 1,mpi_comm_world, istatus, ierr)
731 if (.true. == istrunk)
then
733 else if(.true. == isleaf)
then
736 call errtrp(
'recsolve: never come here')
738 call mpi_send(b, m_pds_ndeg*m_pds_neqns, mpi_real8, imp, 1,mpi_comm_world, ierr)
743end subroutine recsolve
747subroutine trunksolve(b)
764real(8),
dimension(:),
intent(inout) :: b
767real(8),
allocatable,
dimension(:) :: b1, x1, v1, w1, xab1, vtmp1
768real(8),
allocatable,
dimension(:) :: b2, x2, v2, w2, xab2, vtmp2
769real(8),
allocatable,
dimension(:) :: bd, xd
772integer :: neqns, ndeg, neqns1, neqns2, neqnsd, nd
776integer :: i,j,k,l,m,n
779integer,
dimension(MPI_STATUS_SIZE) :: istatus
784if (.false.==m_pds_isready_dc)
then
785 call errtrp(
'trunksolve is not ready.')
799allocate(b1(neqns1*ndeg), x1(neqns1*ndeg), v1(neqns1*ndeg), w1(neqns1*ndeg), xab1(neqns1*ndeg), vtmp1(neqns1*ndeg))
800allocate(b2(neqns2*ndeg), x2(neqns2*ndeg), v2(neqns2*ndeg), w2(neqns2*ndeg), xab2(neqns2*ndeg), vtmp2(neqns2*ndeg))
801allocate(bd(nd), xd(nd))
810call divvec(neqns, ndeg, m_pds_part, b, b1, b2, bd)
815call mpi_send(do_calc, 1, mpi_logical, ip1,1,mpi_comm_world,ierr)
816call mpi_send(do_calc, 1, mpi_logical, ip2,1,mpi_comm_world,ierr)
818call mpi_send(b1, ndeg*neqns1, mpi_real8, ip1,1,mpi_comm_world,ierr)
819call mpi_send(b2, ndeg*neqns2, mpi_real8, ip2,1,mpi_comm_world,ierr)
821call mpi_recv(xab1, ndeg*neqns1, mpi_real8, ip1, 1,mpi_comm_world, istatus, ierr)
822call mpi_recv(xab2, ndeg*neqns2, mpi_real8, ip2, 1,mpi_comm_world, istatus, ierr)
827 bd(k) = bd(k) -(vtmp1(i) * xab1(i))
831 bd(k) = bd(k) -(vtmp2(i) * xab2(i))
842call denssolve(m_pds_d, nd, bd)
855 v1(i) = v1(i) + vtmp1(i)*xd(k)
863 v2(i) = v2(i) + vtmp2(i)*xd(k)
872call mpi_send(do_calc, 1, mpi_logical, ip1,1,mpi_comm_world,ierr)
873call mpi_send(do_calc, 1, mpi_logical, ip2,1,mpi_comm_world,ierr)
875call mpi_send(v1, ndeg*neqns1, mpi_real8, ip1,1,mpi_comm_world,ierr)
876call mpi_send(v2, ndeg*neqns2, mpi_real8, ip2,1,mpi_comm_world,ierr)
878call mpi_recv(w1, ndeg*neqns1, mpi_real8, ip1, 1,mpi_comm_world, istatus, ierr)
879call mpi_recv(w2, ndeg*neqns2, mpi_real8, ip2, 1,mpi_comm_world, istatus, ierr)
894call reovec(neqns, ndeg, m_pds_part, b, x1, x2, xd)
898end subroutine trunksolve
903subroutine leafsolve(b)
907real(8),
dimension(:),
intent(inout) :: b
912if (.false.==m_pds_isready_v)
then
913 call errtrp(
'leafsolve is not ready.')
916call nusol0(b, m_pds_v, ierr)
917if (ierr .ne. 0)
call errtrp(
'leafsolve: nusol0')
920end subroutine leafsolve
935real(8),
allocatable,
dimension(:,:) :: pd
938real(8),
allocatable,
dimension(:) :: vtmp, vpc
939integer :: ndeg, nndeg
940integer :: jcol, iofset, idx
944integer,
dimension(MPI_STATUS_SIZE) :: istatus
948real(8) :: t10s, t10e, t20s, t20e, t30s, t30e
953call elapout(
'calcCtAC: waiting c matrix via MPI')
954call mpi_recv(pc%neqns, 1, mpi_integer, imp, 1,mpi_comm_world, istatus, ierr)
955call mpi_recv(pc%ncol, 1, mpi_integer, imp, 1,mpi_comm_world, istatus, ierr)
956call mpi_recv(pc%nttbr, 1, mpi_integer, imp, 1,mpi_comm_world, istatus, ierr)
957call mpi_recv(pc%ndeg, 1, mpi_integer, imp, 1,mpi_comm_world, istatus, ierr)
963call elapout(
'calcCtAC: begin receive C matrix')
965call mpi_recv(pc%ia, pc%ncol+1, mpi_integer, imp, 1,mpi_comm_world, istatus, ierr)
966call mpi_recv(pc%ja, pc%nttbr, mpi_integer, imp, 1,mpi_comm_world, istatus, ierr)
967call mpi_recv(pc%val, pc%nttbr*nndeg, mpi_real8, imp, 1,mpi_comm_world, istatus, ierr)
969call elapout(
'calcCtAC: end receive C matrix')
971allocate(vtmp(pc%neqns*pc%ndeg))
972allocate(vpc(pc%neqns*pc%ndeg))
977call elapout(
'calcCtAC: start solve CtAC')
985 if (.true. == istrunk)
then
986 call trunksolve(vtmp)
987 else if(.true. == isleaf)
then
990 call errtrp(
'calcCtAC: never come here')
996 jcol = (j+ndeg-1) / ndeg
997 iofset = mod(j+ndeg-1, ndeg)
998 do k=pc%ia(jcol),pc%ia(jcol+1)-1
1001 pd(j,i)=pd(j,i)+pc%val(ndeg*iofset + l,k) * vtmp(ndeg*(idx-1)+l)
1010call elapout(
'calcCtAC: end solve CtAC')
1012call elapout(
'calcCtAC: start send CtAC')
1014call mpi_send(pd, npd*npd, mpi_real8, imp, 1,mpi_comm_world, ierr)
1016deallocate(vtmp,pd,vpc)
1017call elapout(
'calcCtAC: end send CtAC')
1032subroutine mkpart(neqns, nttbr, ndeg, irow, jcol, neqns1, neqns2, neqnsd, part, iperm)
1041integer,
intent(in) :: neqns, nttbr, ndeg
1042integer,
dimension(:),
intent(in) :: irow
1043integer,
dimension(:),
intent(in) :: jcol
1045integer,
intent(out) :: neqns1, neqns2
1046integer,
intent(out) :: neqnsd
1047integer,
dimension(:),
intent(out) :: part
1048integer,
dimension(:),
intent(out) :: iperm
1050integer :: ipart1(neqns),ipart2(neqns),isp(neqns)
1051integer :: loc1, loc2, loc3
1052integer :: i,j,k,l,m,n
1071 if(part(i).eq.1)
then
1075 if(part(i).eq.2)
then
1079 if(part(i).eq.3)
then
1086end subroutine mkpart
1091subroutine divmat(a0, part, iperm, neqns1, neqns2, neqnsd, ndeg, a1, a2, c1, c2, d)
1108integer,
dimension(:),
intent(in) :: part
1109integer,
dimension(:),
intent(in) :: iperm
1110integer,
intent(in) :: neqns1
1111integer,
intent(in) :: neqns2
1112integer,
intent(in) :: neqnsd
1113integer,
intent(in) :: ndeg
1117real(8),
dimension(:,:),
intent(out) :: d
1120integer,
allocatable :: jstat1(:), irpt1(:), irowno1(:)
1121integer,
allocatable :: jstat2(:), irpt2(:), irowno2(:)
1122integer :: nttbr1, nttbr2, ntt1, ntt2
1124integer :: ipass, nndeg
1125integer :: i,j,k,l,m,n
1129allocate(jstat1(neqnsd), jstat2(neqnsd))
1130allocate(irpt1(a0%nttbr), irpt2(a0%nttbr))
1131allocate(irowno1(a0%nttbr), irowno2(a0%nttbr))
1152 if(part(i).eq.1.and.part(j).eq.1)
then
1155 a1%irow(nttbr1)=iperm(i)
1156 a1%jcol(nttbr1)=iperm(j)
1158 a1%val(k,nttbr1)=a0%val(k,l)
1163 if(part(i).eq.2.and.part(j).eq.2)
then
1166 a2%irow(nttbr2)=iperm(i)
1167 a2%jcol(nttbr2)=iperm(j)
1169 a2%val(k,nttbr2)=a0%val(k,l)
1174 if(part(i).eq.1.and.part(j).eq.3)
then
1176 call reserv(iperm(i),iperm(j),jstat1,irpt1,irowno1,ntt1)
1179 call stval(c1,iperm(i),iperm(j),a0%val(1,l),0)
1183 if(part(i).eq.2.and.part(j).eq.3)
then
1185 call reserv(iperm(i),iperm(j),jstat2,irpt2,irowno2,ntt2)
1188 call stval(c2,iperm(i),iperm(j),a0%val(1,l),0)
1192 if(part(i).eq.3.and.part(j).eq.1)
then
1194 call reserv(iperm(j),iperm(i),jstat1,irpt1,irowno1,ntt1)
1197 call stval(c1,iperm(j),iperm(i),a0%val(1,l),1)
1201 if(part(i).eq.3.and.part(j).eq.2)
then
1203 call reserv(iperm(j),iperm(i),jstat2,irpt2,irowno2,ntt2)
1206 call stval(c2,iperm(j),iperm(i),a0%val(1,l),1)
1210 if(part(i).eq.3.and.part(j).eq.3)
then
1214 d(ndeg*(iperm(i)-1)+k,ndeg*(iperm(j)-1)+m)=a0%val(ndeg*(m-1)+k,l)
1216 d(ndeg*(iperm(j)-1)+k,ndeg*(iperm(i)-1)+m)=a0%val(ndeg*(k-1)+m,l)
1223 if(part(i).eq.1.and.part(j).eq.2)
then
1224 write(20,*)
"divmat: wrong",i,j
1227 if(part(i).eq.2.and.part(j).eq.1)
then
1228 write(20,*)
"divmat: wrong",i,j
1242 allocate(a1%irow(nttbr1), a1%jcol(nttbr1), a1%val(ndeg*ndeg,nttbr1))
1243 allocate(a2%irow(nttbr2), a2%jcol(nttbr2), a2%val(ndeg*ndeg,nttbr2))
1249 allocate(c1%ia(c1%ncol+1))
1250 allocate(c1%ja(c1%nttbr))
1251 allocate(c1%val(ndeg*ndeg,c1%nttbr))
1257 allocate(c2%ia(c2%ncol+1))
1258 allocate(c2%ja(c2%nttbr))
1259 allocate(c2%val(ndeg*ndeg,c2%nttbr))
1265 call stiajac(c1,jstat1,irpt1,irowno1)
1266 call stiajac(c2,jstat2,irpt2,irowno2)
1267 deallocate(jstat1,irpt1,irowno1)
1268 deallocate(jstat2,irpt2,irowno2)
1277end subroutine divmat
1282subroutine divvec(neqns, ndeg, part, r, r1, r2, rd)
1286integer,
intent(in) :: neqns, ndeg
1287real(8),
dimension(:),
intent(in) :: r
1288integer,
dimension(:),
intent(in) :: part
1290real(8),
dimension(:),
intent(out) :: r1
1291real(8),
dimension(:),
intent(out) :: r2
1292real(8),
dimension(:),
intent(out) :: rd
1294integer :: idx1, idx2, idxd, idxr
1302 if(part(i).eq.1)
then
1303 r1(idx1:idx1+ndeg-1)=r(idxr:idxr+ndeg-1)
1307 if(part(i).eq.2)
then
1308 r2(idx2:idx2+ndeg-1)=r(idxr:idxr+ndeg-1)
1312 if(part(i).eq.3)
then
1313 rd(idxd:idxd+ndeg-1)=r(idxr:idxr+ndeg-1)
1319end subroutine divvec
1324subroutine reovec(neqns, ndeg, part, r, r1, r2, rd)
1328integer,
intent(in) :: neqns, ndeg
1329integer,
dimension(:),
intent(in) :: part
1330real(8),
dimension(:),
intent(in) :: r1
1331real(8),
dimension(:),
intent(in) :: r2
1332real(8),
dimension(:),
intent(in) :: rd
1334real(8),
dimension(:),
intent(out) :: r
1336integer :: idx1, idx2, idxd, idxr
1345 if(part(i).eq.1)
then
1346 r(idxr:idxr+ndeg-1)=r1(idx1:idx1+ndeg-1)
1350 if(part(i).eq.2)
then
1351 r(idxr:idxr+ndeg-1)=r2(idx2:idx2+ndeg-1)
1355 if(part(i).eq.3)
then
1356 r(idxr:idxr+ndeg-1)=rd(idxd:idxd+ndeg-1)
1362end subroutine reovec
1366subroutine densldu(a,n)
1391integer :: i,j,k,l,m,n
1392real(8),
dimension(n,n) :: a
1403 a(1,j) = a(1,j)/a(1,1)
1408 a(i,i) = a(i,i) - a(i,k)*a(k,i)
1413 a(j,i)=a(j,i)-a(j,k)*a(k,i)
1417 a(i,j)=a(i,j)-a(i,k)*a(k,j)
1419 a(i,j) = a(i,j)/a(i,i)
1426 a(i,j) = a(i,j)/a(j,j)
1438end subroutine densldu
1443subroutine denssolve(a,ndim,r)
1448real(8),
dimension(:,:) :: a
1449real(8),
dimension(:) :: r
1454integer :: i,j,k,l,m,n
1459 r(i)=r(i)-a(i,j)*r(j)
1471 r(i)=r(i)-a(i,j)*r(j)
1475end subroutine denssolve
1491call mpi_comm_size(mpi_comm_world, npe, ierr)
1492call mpi_comm_rank(mpi_comm_world, myid, ierr)
1494if ((mod(npe,2)==0) .and. (myid == npe-1))
then
1510if (myid*2+1 < npe-1)
then
1524subroutine errtrp(mes)
1526write(6,*)
'Error in : process ', myid
1531end subroutine errtrp
subroutine verif0(neqns, ndeg, nttbr, irow, jcol, val, rhs, x)
void version(char *arg)
show version and revision of FrontISTR
void bi_part_directive(int *neqns, int *nttbr, int *irow, int *jcol, int *num_graph1, int *num_graph2, int *num_separator)
subroutine cclsallocation(c)
subroutine m_cclsmatrix_getvec(c, k, v)
subroutine stval(c, i, j, val, itrans)
subroutine reserv(i, j, jstat, irpt, irowno, k)
subroutine stiajac(c, jstat, irpt, irowno)
subroutine, public elapout(mes)
subroutine, public sp_lineq(sp_mat)
void get_part_result(int *num_graph1, int *igraph1, int *num_graph2, int *igraph2, int *num_separator, int *iseparator)