1! This software is Copyright (c) 2012-2015 The Regents of the 2! University of California. All Rights Reserved. 3! 4! Permission to copy and modify this software and its documentation 5! for educational, research and non-profit purposes, without fee, 6! and without a written agreement is hereby granted, provided that 7! the above copyright notice, this paragraph and the following three 8! paragraphs appear in all copies. 9! 10! Permission to make commercial use of this software may be obtained 11! by contacting: 12! 13! Technology Transfer Office 14! 9500 Gilman Drive, Mail Code 0910 15! University of California 16! La Jolla, CA 92093-0910 17! (858) 534-5815 18! invent@ucsd.edu 19! 20! This software program and documentation are copyrighted by The 21! Regents of the University of California. The software program and 22! documentation are supplied "as is", without any accompanying 23! services from The Regents. The Regents does not warrant that the 24! operation of the program will be uninterrupted or error-free. The 25! end-user understands that the program was developed for research 26! purposes and is advised not to rely exclusively on the program for 27! any reason. 28! 29! IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY 30! PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL 31! DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS 32! SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF 33! CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 34! THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY 35! WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 36! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE 37! SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE 38! UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE 39! MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. 40 41 SUBROUTINE NSSOLVER(nFaces, gnNo, dof, nNo, nnz, mynNo, commu, 42 2 cS, face, ls, rowPtr, colPtr, Val, Ri) 43 44 INCLUDE "svLS_STD.h" 45 46 47 INTEGER, INTENT(IN) :: nFaces, gnNo, dof, nNo, nnz, mynNo 48 TYPE(svLS_commuType), INTENT(IN) :: commu 49 TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks) 50 TYPE(svLS_faceType), INTENT(INOUT) :: face(nFaces) 51 TYPE(svLS_lsType), INTENT(INOUT) :: ls 52 INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz) 53 REAL*8, INTENT(IN) :: Val(dof*dof,nnz) 54 REAL*8, INTENT(INOUT) :: Ri(dof,nNo) 55 56 INTEGER i, j, k, iB, iBB, nB, nsd 57 REAL*8 CPUT, NORMS, NORMV, DOTS, DOTV, eps 58 REAL*8, ALLOCATABLE :: U(:,:,:), P(:,:), 59 2 MU(:,:,:), MP(:,:), A(:,:), B(:), xB(:), mK(:,:), mG(:,:), 60 3 mD(:,:), mL(:), Gt(:,:), Rm(:,:), Rc(:), Rmi(:,:), Rci(:) 61 62 nsd = dof - 1 63 iB = ls%RI%mItr 64 nB = 2*iB 65 ALLOCATE(Rm(nsd,nNo), Rc(nNo), Rmi(nsd,nNo), Rci(nNo), 66 2 U(nsd,nNo,iB), P(nNo,iB), MU(nsd,nNo,nB), MP(nNo,nB), 67 3 A(nB,nB), B(nB), xB(nB)) 68 69 Rmi = Ri(1:nsd,:) 70 Rci = Ri(dof,:) 71 72 xB = 0D0 73 B = 0D0 74 Rm = Rmi 75 Rc = Rci 76 eps = SQRT(NORMV(nsd, mynNo, commu, Rm)**2D0 77 2 + NORMS( mynNo, commu, Rc)**2D0) 78 ls%RI%iNorm = eps 79 ls%RI%fNorm = eps 80 ls%CG%callD = 0D0 81 ls%GM%callD = 0D0 82 ls%CG%itr = 0 83 ls%GM%itr = 0 84 ls%RI%callD = CPUT() 85 ls%RI%suc = .FALSE. 86 eps = MAX(ls%RI%absTol,ls%RI%relTol*eps) 87 88 CALL DEPART 89 CALL BCPRE 90 91 DO i=1, ls%RI%mItr 92 iB = 2*i - 1 93 iBB = 2*i 94 ls%RI%dB = ls%RI%fNorm 95 CALL GMRES(nFaces, nsd, nNo, nnz, mynNo, commu, cS, face, 96 2 ls%GM, rowPtr, colPtr, mK, Rm, U(:,:,i)) 97 98 CALL SPARMULVS(nsd, nNo, nnz, commu, cS, 99 2 rowPtr, colPtr, mD, U(:,:,i), P(:,i)) 100 101 P(:,i) = Rc - P(:,i) 102 CALL CGRAD(nFaces, nsd, nNo, nnz, mynNo, commu, cS, face, 103 2 ls%CG, rowPtr, colPtr, Gt, mG, mL, P(:,i)) 104 105 CALL SPARMULSV(nsd, nNo, nnz, commu, cS, 106 2 rowPtr, colPtr, mG, P(:,i), MU(:,:,iB)) 107 108 MU(:,:,iBB) = Rm - MU(:,:,iB) 109 CALL GMRES(nFaces, nsd, nNo, nnz, mynNo, commu, cS, face, 110 2 ls%GM, rowPtr, colPtr, mK, MU(:,:,iBB), U(:,:,i)) 111 112 CALL SPARMULVV(nsd, nNo, nnz, commu, cS, 113 2 rowPtr, colPtr, mK, U(:,:,i), MU(:,:,iBB)) 114 115 CALL ADDBCMUL(BCOP_TYPE_ADD, nFaces, nsd, nNo, mynNo, commu, 116 2 face, U(:,:,i), MU(:,:,iBB)) 117 118 CALL SPARMULSS(nNo, nnz, commu, cS, 119 2 rowPtr, colPtr, mL, P(:,i), MP(:,iB)) 120 121 CALL SPARMULVS(nsd, nNo, nnz, commu, cS, 122 2 rowPtr, colPtr, mD, U(:,:,i), MP(:,iBB)) 123 124 DO k=iB, iBB 125 DO j=1, k - 1 126 A(j,k) = DOTV(nsd, mynNo, commu, MU(:,:,j), MU(:,:,k)) 127 2 + DOTS( mynNo, commu, MP(:,j), MP(:,k)) 128 A(k,j) = A(j,k) 129 END DO 130 A(k,k) = NORMV(nsd, mynNo, commu, MU(:,:,k))**2D0 131 2 + NORMS( mynNo, commu, MP(:,k))**2D0 132 B(k) = DOTV (nsd, mynNo, commu, MU(:,:,k), Rmi) 133 2 + DOTS ( mynNo, commu, MP(:,k), Rci) 134 END DO 135 136 xB = B 137 CALL GE(iBB, A(1:iBB,1:iBB), xB(1:iBB)) 138 139 ls%RI%fNorm = SQRT(ls%RI%iNorm**2D0 - SUM(xB(1:iBB)*B(1:iBB))) 140 IF(ls%RI%fNorm .LT. eps) THEN 141 ls%RI%suc = .TRUE. 142 EXIT 143 END IF 144 145 Rm = Rmi - xB(1)*MU(:,:,1) 146 Rc = Rci - xB(1)*MP(:,1) 147 DO j=2, iBB 148 Rm = Rm - xB(j)*MU(:,:,j) 149 Rc = Rc - xB(j)*MP(:,j) 150 END DO 151 END DO 152 IF (i .GT. ls%RI%mItr) THEN 153 ls%RI%itr = ls%RI%mItr 154 ELSE 155 ls%RI%itr = i 156 157 Rc = Rci - xB(1)*MP(:,1) 158 DO j=2, iBB 159 Rc = Rc - xB(j)*MP(:,j) 160 END DO 161 END IF 162 ls%Resc = NINT(1D2*(NORMS(mynNo, commu, Rc)/ 163 2 ls%RI%fNorm)**2D0) 164 ls%Resm = 100 - ls%Resc 165 166 Rmi = xB(2)*U(:,:,1) 167 Rci = xB(1)*P(:,1) 168 DO i=2, ls%RI%itr 169 iB = 2*i - 1 170 iBB = 2*i 171 172 Rmi = Rmi + xB(iBB)*U(:,:,i) 173 Rci = Rci + xB(iB)*P(:,i) 174 END DO 175 176 ls%RI%callD = CPUT() - ls%RI%callD 177 ls%RI%dB = 1D1*LOG(ls%RI%fNorm/ls%RI%dB) 178 179 Ri(1:nsd,:) = Rmi 180 Ri(dof,:) = Rci 181 182 DEALLOCATE (Rm, Rc, Rmi, Rci, U, P, MU, MP, A, B, mK, mD, mG, mL, 183 2 Gt) 184 185 IF (commu%masF) CALL LOGFILE 186 187 RETURN 188 CONTAINS 189 190!==================================================================== 191 192 SUBROUTINE DEPART 193 194 IMPLICIT NONE 195 196 INTEGER i, j, k, l 197 REAL*8, ALLOCATABLE :: tmp(:) 198 199 ALLOCATE(mK(nsd*nsd,nnz), mG(nsd,nnz), mD(nsd,nnz), mL(nnz), 200 2 Gt(nsd,nnz), tmp((nsd+1)*(nsd+1))) 201 202 IF (nsd .EQ. 2) THEN 203 DO i=1, nnz 204 tmp = Val(:,i) 205 206 mK(1,i) = tmp(1) 207 mK(2,i) = tmp(2) 208 mK(3,i) = tmp(4) 209 mK(4,i) = tmp(5) 210 211 mG(1,i) = tmp(3) 212 mG(2,i) = tmp(6) 213 214 mD(1,i) = tmp(7) 215 mD(2,i) = tmp(8) 216 217 mL(i) = tmp(9) 218 END DO 219 ELSE IF(nsd .EQ. 3) THEN 220 DO i=1, nnz 221 tmp = Val(:,i) 222 223 mK(1,i) = tmp(1) 224 mK(2,i) = tmp(2) 225 mK(3,i) = tmp(3) 226 mK(4,i) = tmp(5) 227 mK(5,i) = tmp(6) 228 mK(6,i) = tmp(7) 229 mK(7,i) = tmp(9) 230 mK(8,i) = tmp(10) 231 mK(9,i) = tmp(11) 232 233 mG(1,i) = tmp(4) 234 mG(2,i) = tmp(8) 235 mG(3,i) = tmp(12) 236 237 mD(1,i) = tmp(13) 238 mD(2,i) = tmp(14) 239 mD(3,i) = tmp(15) 240 241 mL(i) = tmp(16) 242 END DO 243 ELSE 244 PRINT *, "Not defined nsd for DEPART", nsd 245 END IF 246 247 DO i=1, nNo 248 Do j=rowPtr(1,i), rowPtr(2,i) 249 k = colPtr(j) 250 DO l=rowPtr(1,k), rowPtr(2,k) 251 IF (colPtr(l) .EQ. i) THEN 252 Gt(:,l) = -mG(:,j) 253 EXIT 254 END IF 255 END DO 256 END DO 257 END DO 258 259 RETURN 260 END SUBROUTINE DEPART 261 262!==================================================================== 263 264 SUBROUTINE BCPRE 265 266 IMPLICIT NONE 267 268 INTEGER faIn, i, a, Ac 269 REAL*8 NORMV 270 REAL*8, ALLOCATABLE :: v(:,:) 271 272 DO faIn=1, nFaces 273 IF (face(faIn)%coupledFlag) THEN 274 IF (face(faIn)%sharedFlag) THEN 275 IF (.NOT.ALLOCATED(v)) ALLOCATE(v(nsd,nNo)) 276 v = 0D0 277 DO a=1, face(faIn)%nNo 278 Ac = face(faIn)%glob(a) 279 DO i=1, nsd 280 v(i,Ac) = face(faIn)%valM(i,a) 281 END DO 282 END DO 283 face(faIn)%nS = NORMV(nsd, mynNo, commu, v)**2D0 284 ELSE 285 face(faIn)%nS = 0D0 286 DO a=1, face(faIn)%nNo 287 Ac = face(faIn)%glob(a) 288 DO i=1, nsd 289 face(faIn)%nS = face(faIn)%nS + 290 2 face(faIn)%valM(i,a)**2D0 291 END DO 292 END DO 293 END IF 294 END IF 295 END DO 296 297 RETURN 298 END SUBROUTINE BCPRE 299 300!==================================================================== 301 302 SUBROUTINE LOGFILE 303 304 IMPLICIT NONE 305 306 LOGICAL flag 307 INTEGER fid, i, j 308 CHARACTER*16, PARAMETER :: fName = 'svLS_NS.log' 309 310 INQUIRE(FILE=fName, EXIST=flag) 311 312 fid = 11232 313 OPEN(fid, FILE=fName, POSITION='APPEND') 314 315 IF (.NOT.flag) THEN 316 i = 0 317 DO j=1, nFaces 318 IF (face(j)%coupledFlag) i = i + 1 319 END DO 320 WRITE(fid,*) gnNo, commu%nTasks, i 321 END IF 322 323 i = 0 324 IF (ls%RI%suc) i = i + 100 325 IF (ls%GM%suc) i = i + 10 326 IF (ls%CG%suc) i = i + 1 327 328 WRITE(fid,"(I4.3,I3,I4,I5,3I4,3ES9.2E2,3I4)") 329 2 i, ls%RI%itr, ls%GM%itr, ls%CG%itr, 330 3 NINT((ls%RI%CallD-ls%GM%CallD-ls%CG%CallD)/ls%RI%CallD*1D2), 331 4 NINT(ls%GM%callD/ls%RI%CallD*1D2), 332 5 NINT(ls%CG%callD/ls%RI%CallD*1D2), 333 6 ls%RI%iNorm, ls%RI%fNorm/ls%RI%iNorm, ls%RI%CallD, 334 7 ls%Resm, ls%Resc, NINT(ls%RI%dB) 335 336 CLOSE(fid) 337 338 RETURN 339 END SUBROUTINE LOGFILE 340 341 END SUBROUTINE NSSOLVER 342 343