52 VPUBLIC
void Vdpbsl(
double *abd,
int *lda,
int *n,
int *m,
double *b) {
55 int k, kb, la, lb, lm;
59 for (k=1; k<=*n; k++) {
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);
68 for (kb=1; kb<=*n; kb++) {
74 VAT(b, k) /= VAT2(abd, *m+1, k);
76 Vdaxpy(lm, t, RAT2(abd, la, k), 1, RAT(b, lb), 1);
80 VPUBLIC
void Vdaxpy(
int n,
double da,
82 double *dy,
int incy) {
84 int i, ix, iy, m, mp1;
92 if (incx == 1 && incy == 1) {
98 VAT(dy, i) += da * VAT(dx, i);
106 for (i=mp1; i<=n; i+=4) {
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);
117 ix = (-n + 1) * incx + 1;
121 iy = (-n + 1) * incy + 1;
123 for (i=1; i<=n; i++) {
125 VAT(dy, iy) += da * VAT(dx, ix);
133 VPUBLIC
double Vddot(
int n,
double *dx,
int incx,
double *dy,
int incy) {
136 int i, ix, iy, m, mp1;
144 if (incx == 1 && incy == 1) {
151 dtemp += VAT(dx, i) * VAT(dy, i);
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);
171 ix = (-n + 1) * incx + 1;
175 iy = (-n + 1) * incy + 1;
177 for (i=1; i<=n; i++) {
189 VPUBLIC
void Vdpbfa(
double *abd,
int *lda,
int *n,
int *m,
int *info) {
192 int ik, j, jk, k, mu;
198 for(j = 1; j <= *n; j++) {
202 jk = VMAX2(j - *m, 1);
203 mu = VMAX2(*m + 2 - j, 1);
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);
219 s = VAT2(abd, *m + 1, j) - s;
226 VAT2(abd, *m + 1, j) = VSQRT(s);
VPUBLIC void Vdpbsl(double *abd, int *lda, int *n, int *m, double *b)
LINPACK interface.