57VPUBLIC
void Vdpbsl(
double *abd,
int *lda,
int *n,
int *m,
double *b) {
60 int k, kb, la, lb, lm;
64 for (k=1; k<=*n; k++) {
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);
73 for (kb=1; kb<=*n; kb++) {
79 VAT(b, k) /= VAT2(abd, *m+1, k);
81 Vdaxpy(lm, t, RAT2(abd, la, k), 1, RAT(b, lb), 1);
85VPUBLIC
void Vdaxpy(
int n,
double da,
87 double *dy,
int incy) {
89 int i, ix, iy, m, mp1;
97 if (incx == 1 && incy == 1) {
103 VAT(dy, i) += da * VAT(dx, i);
111 for (i=mp1; i<=n; i+=4) {
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);
122 ix = (-n + 1) * incx + 1;
126 iy = (-n + 1) * incy + 1;
128 for (i=1; i<=n; i++) {
130 VAT(dy, iy) += da * VAT(dx, ix);
138VPUBLIC
double Vddot(
int n,
double *dx,
int incx,
double *dy,
int incy) {
141 int i, ix, iy, m, mp1;
149 if (incx == 1 && incy == 1) {
156 dtemp += VAT(dx, i) * VAT(dy, i);
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);
176 ix = (-n + 1) * incx + 1;
180 iy = (-n + 1) * incy + 1;
182 for (i=1; i<=n; i++) {
194VPUBLIC
void Vdpbfa(
double *abd,
int *lda,
int *n,
int *m,
int *info) {
197 int ik, j, jk, k, mu;
203 for(j = 1; j <= *n; j++) {
207 jk = VMAX2(j - *m, 1);
208 mu = VMAX2(*m + 2 - j, 1);
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);
224 s = VAT2(abd, *m + 1, j) - s;
231 VAT2(abd, *m + 1, j) = VSQRT(s);
VPUBLIC void Vdpbsl(double *abd, int *lda, int *n, int *m, double *b)
LINPACK interface.