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 GMRES(nFaces, dof, nNo, nnz, mynNo, commu, cS, face, 42 2 ls, rowPtr, colPtr, Val, R, X) 43 44 INCLUDE "svLS_STD.h" 45 46 INTEGER, INTENT(IN) :: nFaces, dof, nNo, nnz, mynNo 47 TYPE(svLS_commuType), INTENT(IN) :: commu 48 TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks) 49 TYPE(svLS_faceType), INTENT(IN) :: face(nFaces) 50 TYPE(svLS_subLsType), INTENT(INOUT) :: ls 51 INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz) 52 REAL*8, INTENT(IN) :: Val(dof*dof,nnz), R(dof,nNo) 53 REAL*8, INTENT(OUT) :: X(dof,nNo) 54 55 INTEGER i, j, k, l 56 REAL*8 CPUT, NORMV, DOTV 57 REAL*8 eps, tmp, time 58 REAL*8, ALLOCATABLE :: u(:,:,:), h(:,:), unCondU(:,:), y(:), c(:), 59 2 s(:), err(:) 60 61 ALLOCATE(h(ls%sD+1,ls%sD), u(dof,nNo,ls%sD+1), unCondU(dof,nNo), 62 2 y(ls%sD), c(ls%sD), s(ls%sD), err(ls%sD+1)) 63 64 time = CPUT() 65 ls%suc = .FALSE. 66 67 X = 0D0 68 DO l=1, ls%mItr 69 IF (l .EQ. 1) THEN 70 u(:,:,1) = R 71 ELSE 72 ls%itr = ls%itr + 1 73 CALL SPARMULVV(dof, nNo, nnz, commu, cS, rowPtr, colPtr, 74 2 Val, X, u(:,:,1)) 75 76 CALL ADDBCMUL(BCOP_TYPE_ADD, nFaces, dof, nNo, mynNo, commu, 77 2 face, X, u(:,:,1)) 78 79 u(:,:,1) = R - u(:,:,1) 80 END IF 81 IF (ANY(face%coupledFlag)) THEN 82 unCondU = u(:,:,1) 83 CALL ADDBCMUL(BCOP_TYPE_PRE, nFaces, dof, nNo, mynNo, commu, 84 2 face, unCondU, u(:,:,1)) 85 END IF 86 87 err(1) = NORMV(dof, mynNo, commu, u(:,:,1)) 88 IF (l .EQ. 1) THEN 89 eps = err(1) 90 IF (eps .LE. ls%absTol) THEN 91 ls%callD = 0D0 92 ls%dB = 0D0 93 RETURN 94 END IF 95 ls%iNorm = eps 96 ls%fNorm = eps 97 eps = MAX(ls%absTol,ls%relTol*eps) 98 END IF 99 ls%dB = ls%fNorm 100 u(:,:,1) = u(:,:,1)/err(1) 101 DO i=1, ls%sD 102 ls%itr = ls%itr + 1 103 CALL SPARMULVV(dof, nNo, nnz, commu, cS, rowPtr, colPtr, 104 2 Val, u(:,:,i), u(:,:,i+1)) 105 106 CALL ADDBCMUL(BCOP_TYPE_ADD, nFaces, dof, nNo, mynNo, commu, 107 2 face, u(:,:,i), u(:,:,i+1)) 108 109 IF (ANY(face%coupledFlag)) THEN 110 unCondU = u(:,:,i+1) 111 CALL ADDBCMUL(BCOP_TYPE_PRE, nFaces, dof, nNo, mynNo, 112 2 commu, face, unCondU, u(:,:,i+1)) 113 END IF 114 DO j=1, i 115 h(j,i) = DOTV(dof, mynNo, commu, u(:,:,i+1), u(:,:,j)) 116 u(:,:,i+1) = u(:,:,i+1) - h(j,i)*u(:,:,j) 117 END DO 118 h(i+1,i) = NORMV(dof, mynNo, commu, u(:,:,i+1)) 119 120 u(:,:,i+1) = u(:,:,i+1)/h(i+1,i) 121 DO j=1, i-1 122 tmp = c(j)*h(j,i) + s(j)*h(j+1,i) 123 h(j+1,i) = -s(j)*h(j,i) + c(j)*h(j+1,i) 124 h(j,i) = tmp 125 END DO 126 tmp = SQRT(h(i,i)*h(i,i) + h(i+1,i)*h(i+1,i)) 127 c(i) = h(i,i)/tmp 128 s(i) = h(i+1,i)/tmp 129 h(i,i) = tmp 130 h(i+1,i) = 0D0 131 err(i+1) = -s(i)*err(i) 132 err(i) = c(i)*err(i) 133 IF (ABS(err(i+1)) .LT. eps) THEN 134 ls%suc = .TRUE. 135 EXIT 136 END IF 137 END DO 138 IF (i .GT. ls%sD) i = ls%sD 139 140 y = err(1:i) 141 DO j=i, 1, -1 142 DO k=j+1, i 143 y(j) = y(j) - h(j,k)*y(k) 144 END DO 145 y(j) = y(j)/h(j,j) 146 END DO 147 148 DO j=1, i 149 X = X + u(:,:,j)*y(j) 150 END DO 151 ls%fNorm = ABS(err(i+1)) 152 IF (ls%suc) EXIT 153 END DO 154 155 ls%callD = CPUT() - time + ls%callD 156 ls%dB = 1D1*LOG(ls%fNorm/ls%dB) 157 158 RETURN 159 END SUBROUTINE GMRES 160 161!==================================================================== 162 163 SUBROUTINE GMRESV(nFaces, dof, nNo, nnz, mynNo, commu, cS, face, 164 2 ls, rowPtr, colPtr, Val, R) 165 166 INCLUDE "svLS_STD.h" 167 168 INTEGER, INTENT(IN) :: nFaces, dof, nNo, nnz, mynNo 169 TYPE(svLS_commuType), INTENT(IN) :: commu 170 TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks) 171 TYPE(svLS_faceType), INTENT(IN) :: face(nFaces) 172 TYPE(svLS_subLsType), INTENT(INOUT) :: ls 173 INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz) 174 REAL*8, INTENT(IN) :: Val(dof*dof,nnz) 175 REAL*8, INTENT(INOUT) :: R(dof,nNo) 176 177 INTEGER i, j, k, l 178 REAL*8 CPUT, NORMV, DOTV 179 REAL*8 eps, tmp 180 REAL*8, ALLOCATABLE :: u(:,:,:), h(:,:), X(:,:), y(:), c(:), s(:), 181 2 err(:) 182 183 ALLOCATE(h(ls%sD+1,ls%sD), u(dof,nNo,ls%sD+1), X(dof,nNo), 184 2 y(ls%sD), c(ls%sD), s(ls%sD), err(ls%sD+1)) 185 186 ls%callD = CPUT() 187 ls%suc = .FALSE. 188 eps = NORMV(dof, mynNo, commu, R) 189 ls%iNorm = eps 190 ls%fNorm = eps 191 eps = MAX(ls%absTol,ls%relTol*eps) 192 ls%itr = 0 193 X = 0D0 194 195 IF (ls%iNorm .LE. ls%absTol) THEN 196 ls%callD = 0D0 197 ls%dB = 0D0 198 RETURN 199 END IF 200 DO l=1, ls%mItr 201 ls%dB = ls%fNorm 202 ls%itr = ls%itr + 1 203 CALL SPARMULVV(dof, nNo, nnz, commu, cS, rowPtr, colPtr, 204 2 Val, X, u(:,:,1)) 205 CALL ADDBCMUL(BCOP_TYPE_ADD, nFaces, dof, nNo, mynNo, commu, 206 2 face, X, u(:,:,1)) 207 208 u(:,:,1) = R - u(:,:,1) 209 err(1) = NORMV(dof, mynNo, commu, u(:,:,1)) 210 u(:,:,1) = u(:,:,1)/err(1) 211 DO i=1, ls%sD 212 ls%itr = ls%itr + 1 213 CALL SPARMULVV(dof, nNo, nnz, commu, cS, rowPtr, colPtr, 214 2 Val, u(:,:,i), u(:,:,i+1)) 215 CALL ADDBCMUL(BCOP_TYPE_ADD, nFaces, dof, nNo, mynNo, commu, 216 2 face, u(:,:,i), u(:,:,i+1)) 217 218 DO j=1, i 219 h(j,i) = DOTV(dof, mynNo, commu, u(:,:,i+1), u(:,:,j)) 220 u(:,:,i+1) = u(:,:,i+1) - h(j,i)*u(:,:,j) 221 END DO 222 h(i+1,i) = NORMV(dof, mynNo, commu, u(:,:,i+1)) 223 u(:,:,i+1) = u(:,:,i+1)/h(i+1,i) 224 DO j=1, i-1 225 tmp = c(j)*h(j,i) + s(j)*h(j+1,i) 226 h(j+1,i) = -s(j)*h(j,i) + c(j)*h(j+1,i) 227 h(j,i) = tmp 228 END DO 229 tmp = SQRT(h(i,i)*h(i,i) + h(i+1,i)*h(i+1,i)) 230 c(i) = h(i,i)/tmp 231 s(i) = h(i+1,i)/tmp 232 h(i,i) = tmp 233 h(i+1,i) = 0D0 234 err(i+1) = -s(i)*err(i) 235 err(i) = c(i)*err(i) 236 IF (ABS(err(i+1)) .LT. eps) THEN 237 ls%suc = .TRUE. 238 EXIT 239 END IF 240 END DO 241 IF (i .GT. ls%sD) i = ls%sD 242 243 y = err(1:i) 244 DO j=i, 1, -1 245 DO k=j+1, i 246 y(j) = y(j) - h(j,k)*y(k) 247 END DO 248 y(j) = y(j)/h(j,j) 249 END DO 250 251 DO j=1, i 252 X = X + u(:,:,j)*y(j) 253 END DO 254 ls%fNorm = ABS(err(i+1)) 255 IF (ls%suc) EXIT 256 END DO 257 R = X 258 ls%callD = CPUT() - ls%callD 259 ls%dB = 1D1*LOG(ls%fNorm/ls%dB) 260 261 RETURN 262 END SUBROUTINE GMRESV 263 264