1 /**
2 * @ingroup PMGC
3 * @author Tucker Beck [fortran ->c translation], Michael Holst [original]
4 * @brief
5 * @version $Id:
6 *
7 * @attention
8 * @verbatim
9 *
10 * APBS -- Adaptive Poisson-Boltzmann Solver
11 *
12 * Nathan A. Baker (nathan.baker@pnl.gov)
13 * Pacific Northwest National Laboratory
14 *
15 * Additional contributing authors listed in the code documentation.
16 *
17 * Copyright (c) 2010-2014 Battelle Memorial Institute. Developed at the Pacific Northwest National Laboratory, operated by Battelle Memorial Institute, Pacific Northwest Division for the U.S. Department Energy. Portions Copyright (c) 2002-2010, Washington University in St. Louis. Portions Copyright (c) 2002-2010, Nathan A. Baker. Portions Copyright (c) 1999-2002, The Regents of the University of California. Portions Copyright (c) 1995, Michael Holst.
18 * All rights reserved.
19 *
20 *
21 * Redistribution and use in source and binary forms, with or without
22 * modification, are permitted provided that the following conditions are met:
23 *
24 * - Redistributions of source code must retain the above copyright notice, this
25 * list of conditions and the following disclaimer.
26 *
27 * - Redistributions in binary form must reproduce the above copyright notice,
28 * this list of conditions and the following disclaimer in the documentation
29 * and/or other materials provided with the distribution.
30 *
31 * - Neither the name of Washington University in St. Louis nor the names of its
32 * contributors may be used to endorse or promote products derived from this
33 * software without specific prior written permission.
34 *
35 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
36 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
37 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
38 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
39 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
40 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
41 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
42 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
43 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
44 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
45 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
46 *
47 * @endverbatim
48 */
49
50 #include "mlinpckd.h"
51
Vdpbsl(double * abd,int * lda,int * n,int * m,double * b)52 VPUBLIC void Vdpbsl(double *abd, int *lda, int *n, int *m, double *b) {
53
54 double t;
55 int k, kb, la, lb, lm;
56
57 MAT2(abd, *lda, 1);
58
59 for (k=1; k<=*n; k++) {
60 lm = VMIN2(k-1, *m);
61 la = *m + 1 - lm;
62 lb = k - lm;
63 t = Vddot(lm, RAT2(abd, la, k ), 1, RAT(b, lb), 1 );
64 VAT(b, k) = (VAT(b, k) - t) / VAT2(abd, *m+1, k);
65 }
66
67 // Solve R*X = Y
68 for (kb=1; kb<=*n; kb++) {
69
70 k = *n + 1 - kb;
71 lm = VMIN2(k-1, *m);
72 la = *m + 1 - lm;
73 lb = k - lm;
74 VAT(b, k) /= VAT2(abd, *m+1, k);
75 t = -VAT(b, k);
76 Vdaxpy(lm, t, RAT2(abd, la, k), 1, RAT(b, lb), 1);
77 }
78 }
79
Vdaxpy(int n,double da,double * dx,int incx,double * dy,int incy)80 VPUBLIC void Vdaxpy(int n, double da,
81 double *dx, int incx,
82 double *dy, int incy) {
83
84 int i, ix, iy, m, mp1;
85
86 if (n <= 0)
87 return;
88
89 if (da == 0)
90 return;
91
92 if (incx == 1 && incy == 1) {
93
94 m = n % 4;
95 if (m != 0) {
96
97 for (i=1; i<=m; i++)
98 VAT(dy, i) += da * VAT(dx, i);
99 }
100
101 if (n < 4)
102 return;
103
104 mp1 = m + 1;
105
106 for (i=mp1; i<=n; i+=4) {
107
108 VAT(dy, i ) += da * VAT(dx, i );
109 VAT(dy, i+1) += da * VAT(dx, i+1);
110 VAT(dy, i+2) += da * VAT(dx, i+2);
111 VAT(dy, i+3) += da * VAT(dx, i+3);
112 }
113 } else {
114
115 ix = 1;
116 if (incx < 0 )
117 ix = (-n + 1) * incx + 1;
118
119 iy = 1;
120 if (incy < 0 )
121 iy = (-n + 1) * incy + 1;
122
123 for (i=1; i<=n; i++) {
124
125 VAT(dy, iy) += da * VAT(dx, ix);
126 ix += incx;
127 iy += incy;
128 }
129 }
130 }
131
132
Vddot(int n,double * dx,int incx,double * dy,int incy)133 VPUBLIC double Vddot(int n, double *dx, int incx, double *dy, int incy) {
134
135 double dtemp;
136 int i, ix, iy, m, mp1;
137
138 double ddot = 0.0;
139 dtemp = 0.0;
140
141 if (n <= 0)
142 return ddot;
143
144 if (incx == 1 && incy == 1) {
145
146 m = n % 5;
147
148 if (m != 0) {
149
150 for (i=1; i<=m; i++)
151 dtemp += VAT(dx, i) * VAT(dy, i);
152
153 if (n < 5) {
154 ddot = dtemp;
155 return ddot;
156 }
157 }
158
159 mp1 = m + 1;
160
161 for (i=mp1; i<=n; i+=5)
162 dtemp += VAT(dx, i) * VAT(dy, i)
163 + VAT(dx, i+1) * VAT(dy, i+1)
164 + VAT(dx, i+2) * VAT(dy, i+2)
165 + VAT(dx, i+3) * VAT(dy, i+3)
166 + VAT(dx, i+4) * VAT(dy, i+4);
167 } else {
168
169 ix = 1;
170 if (incx < 0)
171 ix = (-n + 1) * incx + 1;
172
173 iy = 1;
174 if (incy < 0)
175 iy = (-n + 1) * incy + 1;
176
177 for (i=1; i<=n; i++) {
178 ix += incx;
179 iy += incy;
180 }
181 }
182
183 ddot = dtemp;
184 return ddot;
185 }
186
187
188
Vdpbfa(double * abd,int * lda,int * n,int * m,int * info)189 VPUBLIC void Vdpbfa(double *abd, int *lda, int *n, int *m, int *info) {
190
191 double t, s;
192 int ik, j, jk, k, mu;
193
194 MAT2(abd, *lda, 1);
195
196 *info = 0;
197
198 for(j = 1; j <= *n; j++) {
199
200 s = 0.0;
201 ik = *m + 1;
202 jk = VMAX2(j - *m, 1);
203 mu = VMAX2(*m + 2 - j, 1);
204
205 if (*m >= mu ) {
206
207 for(k = mu; k <= *m; k++) {
208 t = VAT2(abd, k, j) - Vddot(k - mu,
209 RAT2(abd, ik, jk), 1,
210 RAT2(abd, mu, j), 1);
211 t /= VAT2(abd, *m + 1, jk);
212 VAT2(abd, k, j) = t;
213 s += t * t;
214 ik--;
215 jk++;
216 }
217 }
218
219 s = VAT2(abd, *m + 1, j) - s;
220
221 if (s <= 0.0) {
222 *info = j;
223 break;
224 }
225
226 VAT2(abd, *m + 1, j) = VSQRT(s);
227 }
228 }
229