xref: /phasta/svLS/CGRAD.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)
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 CGRAD(nFaces, dof, nNo, nnz, mynNo, commu, cS, face,
42     2   ls, rowPtr, colPtr, D, G, L, R)
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) :: D(dof,nnz), G(dof,nnz), L(nnz)
53      REAL*8, INTENT(INOUT) :: R(nNo)
54
55      INTEGER i
56      REAL*8 errO, err, alpha, eps, time
57      REAL*8 CPUT, NORMS, DOTS
58      REAL*8, ALLOCATABLE :: X(:), P(:), SP(:), DGP(:), GP(:,:),
59     2   unCondU(:,:)
60
61      ALLOCATE(X(nNo), P(nNo), SP(nNo), DGP(nNo), GP(dof,nNo),
62     2   unCondU(dof,nNo))
63
64      time     = CPUT()
65      ls%suc   = .FALSE.
66      errO     = NORMS(mynNo, commu, R)
67      ls%iNorm = errO
68      eps      = MAX(ls%absTol,ls%relTol*errO)
69      eps      = eps*eps
70      X        = 0D0
71
72      err      = errO
73      err      = err*err
74      P        = R
75      DO i=1, ls%mItr
76         IF (err .LT. eps) THEN
77            ls%suc = .TRUE.
78            EXIT
79         END IF
80         errO = err
81         CALL SPARMULSV(dof, nNo, nnz, commu, cS, rowPtr, colPtr,
82     2      G, P, GP)
83
84         IF (ANY(face%coupledFlag)) THEN
85            unCondU = GP
86            CALL ADDBCMUL(BCOP_TYPE_PRE, nFaces, dof, nNo, mynNo, commu,
87     2         face, unCondU, GP)
88         END IF
89
90         CALL SPARMULVS(dof, nNo, nnz, commu, cS, rowPtr, colPtr,
91     2      D, GP, DGP)
92
93         CALL SPARMULSS(     nNo, nnz, commu, cS, rowPtr, colPtr,
94     2      L, P, SP)
95
96         SP    = SP - DGP
97         alpha = errO/DOTS(mynNo, commu, P, SP)
98         X     = X + alpha*P
99         R     = R - alpha*SP
100         err   = NORMS(mynNo, commu, R)
101         err   = err*err
102         P     = R + err/errO*P
103      END DO
104      R        = X
105      ls%fNorm = SQRT(err)
106      ls%callD = CPUT() - time + ls%callD
107      ls%dB    = 5D0*LOG(err/errO)
108      ls%itr   = ls%itr + i - 1
109
110      DEALLOCATE(X, P, SP, DGP, GP)
111
112      RETURN
113      END SUBROUTINE CGRAD
114
115!====================================================================
116
117      SUBROUTINE CGRADV(dof, nNo, nnz, mynNo, commu, cS, ls,
118     2   rowPtr, colPtr, K, R)
119
120      INCLUDE "svLS_STD.h"
121
122      INTEGER, INTENT(IN) :: dof, nNo, nnz, mynNo
123      TYPE(svLS_commuType), INTENT(IN) :: commu
124      TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks)
125      TYPE(svLS_subLsType), INTENT(INOUT) :: ls
126      INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz)
127      REAL*8, INTENT(IN) :: K(dof*dof,nnz)
128      REAL*8, INTENT(INOUT) :: R(dof,nNo)
129
130      INTEGER i
131      REAL*8 errO, err, alpha, eps
132      REAL*8 CPUT, NORMV, DOTV
133      REAL*8, ALLOCATABLE :: P(:,:), KP(:,:), X(:,:)
134
135      ALLOCATE(P(dof,nNo), KP(dof,nNo), X(dof,nNo))
136
137      ls%callD = CPUT()
138      ls%suc   = .FALSE.
139      err      = NORMV(dof, mynNo, commu, R)
140      ls%iNorm = err
141      eps      = MAX(ls%absTol,ls%relTol*err)
142      eps      = eps*eps
143      err      = err*err
144      X        = 0D0
145      P        = R
146
147      DO i=1, ls%mItr
148         IF (err .LT. eps) THEN
149            ls%suc = .TRUE.
150            EXIT
151         END IF
152         errO = err
153         CALL SPARMULVV(dof, nNo, nnz, commu, cS, rowPtr, colPtr,
154     2      K, P, KP)
155
156         alpha = errO/DOTV(dof, mynNo, commu, P, KP)
157         X     = X + alpha*P
158         R     = R - alpha*KP
159         err   = NORMV(dof, mynNo, commu, R)
160         err   = err*err
161         P = R + err/errO*P
162      END DO
163
164      R        = X
165      ls%itr   = i - 1
166      ls%fNorm = SQRT(err)
167      ls%callD = CPUT() - ls%callD
168      ls%dB    = 5D0*LOG(err/errO)
169      IF (i .GT. ls%mItr) ls%itr = ls%mItr
170
171      RETURN
172      END SUBROUTINE CGRADV
173
174!====================================================================
175
176      SUBROUTINE CGRADS(nNo, nnz, mynNo, commu, cS, ls,
177     2   rowPtr, colPtr, K, R)
178
179      INCLUDE "svLS_STD.h"
180
181      INTEGER, INTENT(IN) :: nNo, nnz, mynNo
182      TYPE(svLS_commuType), INTENT(IN) :: commu
183      TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks)
184      TYPE(svLS_subLsType), INTENT(INOUT) :: ls
185      INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz)
186      REAL*8, INTENT(IN) :: K(nnz)
187      REAL*8, INTENT(INOUT) :: R(nNo)
188
189      INTEGER i
190      REAL*8 errO, err, alpha, eps
191      REAL*8 CPUT, NORMS, DOTS
192      REAL*8, ALLOCATABLE :: P(:), KP(:), X(:)
193
194      ALLOCATE(P(nNo), KP(nNo), X(nNo))
195
196      ls%callD = CPUT()
197      ls%suc   = .FALSE.
198      err      = NORMS(mynNo, commu, R)
199      ls%iNorm = err
200      eps      = MAX(ls%absTol,ls%relTol*err)
201      eps      = eps*eps
202      err      = err*err
203      X        = 0D0
204      P        = R
205
206      DO i=1, ls%mItr
207         IF (err .LT. eps) THEN
208            ls%suc = .TRUE.
209            EXIT
210         END IF
211         errO = err
212         CALL SPARMULSS(nNo, nnz, commu, cS, rowPtr, colPtr, K, P, KP)
213         alpha = errO/DOTS(mynNo, commu, P, KP)
214         X     = X + alpha*P
215         R     = R - alpha*KP
216         err   = NORMS(mynNo, commu, R)
217         err   = err*err
218         P = R + err/errO*P
219      END DO
220
221      R        = X
222      ls%itr   = i - 1
223      ls%fNorm = SQRT(err)
224      ls%callD = CPUT() - ls%callD
225      ls%dB    = 5D0*LOG(err/errO)
226      IF (i .GT. ls%mItr) ls%itr = ls%mItr
227
228      RETURN
229      END SUBROUTINE CGRADS
230