1 subroutine SolDir (y, ac, yold, acold, 2 & x, 3 & iBC, BC, 4 & res, 5 & iper, ilwork, 6 & shp, shgl, shpb, shglb) 7c 8c---------------------------------------------------------------------- 9c direct solver 10c---------------------------------------------------------------------- 11c 12 use pointer_data 13 14 include "common.h" 15 include "mpif.h" 16 include "auxmpi.h" 17c 18 dimension y(nshg,ndof), ac(nshg,ndof), 19 & yold(nshg,ndof), acold(nshg,ndof), 20 & x(numnp,nsd), 21 & iBC(nshg), BC(nshg,ndofBC), 22 & res(nshg,nflow), 23 & ilwork(nlwork), iper(nshg) 24c 25 dimension shp(MAXTOP,maxsh,MAXQPT), 26 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 27 & shpb(MAXTOP,maxsh,MAXQPT), 28 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 29 30 real*8 yAlpha(nshg,5), acAlpha(nshg,5) 31 32c 33 dimension dyf(4*nshg), indx(4*nshg), solinc(nshg,4) 34 35 dimension globMas(4*nshg,4*nshg) 36 37 write (*,*) 'Warning: using direct solver...' 38c 39c.... set the element matrix flag 40c 41 lhs = 1 ! always 42c 43c.... compute solution at n+alpha 44c 45 call itrYAlpha( yold, acold, y, ac, yAlpha, acAlpha) 46c 47c.... *******************>> Element Data Formation <<****************** 48c 49c.... form the LHS matrices, the residual vector 50c 51 call ElmGMR (yAlpha, acAlpha, x, 52 & shp, shgl, iBC, 53 & BC, shpb, shglb, 54 & res, iper, ilwork, 55 & rowp, colm, lhsK, 56 & lhsP, rerr ) 57 58 globMas = zero 59 npro = numel 60c 61cccc need to assemble here! 62c 63 call bc3Global(globMas, iBC) 64c 65cDEBUG: write the global matrix (nonzero blocks) 66c 67 do i=1,nshg 68 i0 = (i-1)*4 69 70 do j=1,nshg 71 j0 = (j-1)*4 72 73 if (globMas(i0+1,j0+1) .ne. 0) then 74 write (544,21) i,j 75 do ii=1,4 76 write(543,20) (globMas(i0+ii,j0+kk), kk=1,4) 77 enddo 78 endif 79 enddo 80 enddo 81 82 20 format (4(2x,e14.7)) 83 21 format (2(2x,i8)) 84 85c$$$ stop 86 87c 88c.... LU factor the mass matrix 89c 90 indx = 0 91 call ludcmp(globMas, 4*nshg, 4*nshg, indx, d) 92 write(543,*) 'rhs' 93 do i=1, nshg 94 i0 = 4*(i-1) 95 dyf(i0+1) = res(i,1) 96 dyf(i0+2) = res(i,2) 97 dyf(i0+3) = res(i,3) 98 dyf(i0+4) = res(i,4) 99 write(543,20) (dyf(i0+j),j=1,4) 100 enddo 101c 102c.... back-substitute to find dY 103c 104 call lubksb(globMas, 4*nshg, 4*nshg, indx, dyf) 105c 106 write(543,*) 'soln' 107 108 do i=1,nshg 109 i0 = 4*(i-1) 110 solinc(i,1) = dyf(i0+1) 111 solinc(i,2) = dyf(i0+2) 112 solinc(i,3) = dyf(i0+3) 113 solinc(i,4) = dyf(i0+4) 114 enddo 115c 116c 117c.... Now, you satisfy the boundary conditions to newly 118c obtained p,u,v,w 119c 120c 121c You have to set boundary conditions first so Dy distributes 122c 123 call itrCorrect ( y, ac, solinc, iBC) 124 call itrBC (y, ac, iBC, BC, iper, ilwork) 125c 126c.... output the statistics 127c 128 call rstatic (res, y, solinc) 129c 130c.... end 131c 132 return 133 end 134