APBS  3.0.0
mlinpckd.c
1 
55 #include "mlinpckd.h"
56 
57 VPUBLIC 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 
85 VPUBLIC 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 
138 VPUBLIC 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 
194 VPUBLIC 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