APBS 3.0.0
Loading...
Searching...
No Matches
mlinpckd.c
1
55#include "mlinpckd.h"
56
57VPUBLIC void Vdpbsl(double *abd, int *lda, int *n, int *m, double *b) {
58
59 double t;
60 int k, kb, la, lb, lm;
61
62 MAT2(abd, *lda, 1);
63
64 for (k=1; k<=*n; k++) {
65 lm = VMIN2(k-1, *m);
66 la = *m + 1 - lm;
67 lb = k - lm;
68 t = Vddot(lm, RAT2(abd, la, k ), 1, RAT(b, lb), 1 );
69 VAT(b, k) = (VAT(b, k) - t) / VAT2(abd, *m+1, k);
70 }
71
72 // Solve R*X = Y
73 for (kb=1; kb<=*n; kb++) {
74
75 k = *n + 1 - kb;
76 lm = VMIN2(k-1, *m);
77 la = *m + 1 - lm;
78 lb = k - lm;
79 VAT(b, k) /= VAT2(abd, *m+1, k);
80 t = -VAT(b, k);
81 Vdaxpy(lm, t, RAT2(abd, la, k), 1, RAT(b, lb), 1);
82 }
83}
84
85VPUBLIC void Vdaxpy(int n, double da,
86 double *dx, int incx,
87 double *dy, int incy) {
88
89 int i, ix, iy, m, mp1;
90
91 if (n <= 0)
92 return;
93
94 if (da == 0)
95 return;
96
97 if (incx == 1 && incy == 1) {
98
99 m = n % 4;
100 if (m != 0) {
101
102 for (i=1; i<=m; i++)
103 VAT(dy, i) += da * VAT(dx, i);
104 }
105
106 if (n < 4)
107 return;
108
109 mp1 = m + 1;
110
111 for (i=mp1; i<=n; i+=4) {
112
113 VAT(dy, i ) += da * VAT(dx, i );
114 VAT(dy, i+1) += da * VAT(dx, i+1);
115 VAT(dy, i+2) += da * VAT(dx, i+2);
116 VAT(dy, i+3) += da * VAT(dx, i+3);
117 }
118 } else {
119
120 ix = 1;
121 if (incx < 0 )
122 ix = (-n + 1) * incx + 1;
123
124 iy = 1;
125 if (incy < 0 )
126 iy = (-n + 1) * incy + 1;
127
128 for (i=1; i<=n; i++) {
129
130 VAT(dy, iy) += da * VAT(dx, ix);
131 ix += incx;
132 iy += incy;
133 }
134 }
135}
136
137
138VPUBLIC double Vddot(int n, double *dx, int incx, double *dy, int incy) {
139
140 double dtemp;
141 int i, ix, iy, m, mp1;
142
143 double ddot = 0.0;
144 dtemp = 0.0;
145
146 if (n <= 0)
147 return ddot;
148
149 if (incx == 1 && incy == 1) {
150
151 m = n % 5;
152
153 if (m != 0) {
154
155 for (i=1; i<=m; i++)
156 dtemp += VAT(dx, i) * VAT(dy, i);
157
158 if (n < 5) {
159 ddot = dtemp;
160 return ddot;
161 }
162 }
163
164 mp1 = m + 1;
165
166 for (i=mp1; i<=n; i+=5)
167 dtemp += VAT(dx, i) * VAT(dy, i)
168 + VAT(dx, i+1) * VAT(dy, i+1)
169 + VAT(dx, i+2) * VAT(dy, i+2)
170 + VAT(dx, i+3) * VAT(dy, i+3)
171 + VAT(dx, i+4) * VAT(dy, i+4);
172 } else {
173
174 ix = 1;
175 if (incx < 0)
176 ix = (-n + 1) * incx + 1;
177
178 iy = 1;
179 if (incy < 0)
180 iy = (-n + 1) * incy + 1;
181
182 for (i=1; i<=n; i++) {
183 ix += incx;
184 iy += incy;
185 }
186 }
187
188 ddot = dtemp;
189 return ddot;
190}
191
192
193
194VPUBLIC void Vdpbfa(double *abd, int *lda, int *n, int *m, int *info) {
195
196 double t, s;
197 int ik, j, jk, k, mu;
198
199 MAT2(abd, *lda, 1);
200
201 *info = 0;
202
203 for(j = 1; j <= *n; j++) {
204
205 s = 0.0;
206 ik = *m + 1;
207 jk = VMAX2(j - *m, 1);
208 mu = VMAX2(*m + 2 - j, 1);
209
210 if (*m >= mu ) {
211
212 for(k = mu; k <= *m; k++) {
213 t = VAT2(abd, k, j) - Vddot(k - mu,
214 RAT2(abd, ik, jk), 1,
215 RAT2(abd, mu, j), 1);
216 t /= VAT2(abd, *m + 1, jk);
217 VAT2(abd, k, j) = t;
218 s += t * t;
219 ik--;
220 jk++;
221 }
222 }
223
224 s = VAT2(abd, *m + 1, j) - s;
225
226 if (s <= 0.0) {
227 *info = j;
228 break;
229 }
230
231 VAT2(abd, *m + 1, j) = VSQRT(s);
232 }
233}
VPUBLIC void Vdpbsl(double *abd, int *lda, int *n, int *m, double *b)
LINPACK interface.
Definition mlinpckd.c:57