APBS  1.5
mlinpckd.c
1 
50 #include "mlinpckd.h"
51 
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 
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 
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 
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 }
VPUBLIC void Vdpbsl(double *abd, int *lda, int *n, int *m, double *b)
LINPACK interface.
Definition: mlinpckd.c:52