Line data Source code
1 : #include <stdio.h>
2 : #include <tgmath.h>
3 :
4 : #include "dist.h"
5 : #include "mathext.h"
6 : #include "qrfactorization.h"
7 :
8 : #define map(i, j) ((j)-1) * lda + ((i)-1)
9 : #define a(i, j) a[ map(i, j) ]
10 : #define q(i, j) q[ map(i, j) ]
11 : #define map1(i) (i) - 1
12 : #define x(i) x[ map1(i) ]
13 : #define y(i) y[ map1(i) ]
14 : #define ipiv(i) ipiv[ map1(i) ]
15 : #define wrk(i) wrk[ map1(i) ]
16 : #define tau(i) tau[ map1(i) ]
17 : #define swap(T, x, y) \
18 : { \
19 : T tmp = (x); \
20 : x = (y); \
21 : y = tmp; \
22 : }
23 :
24 : // Compute the QR factorization with column pivoting on all columns
25 : // of an M by N matrix A
26 : // A*P = Q*R
27 : //
28 : // The matrix Q is represented as a product of elementary reflectors
29 : //
30 : // Q = H(1) H(2) . . . H(k), where k = min(m,n).
31 : //
32 : // Each H(i) has the form
33 : //
34 : // H(i) = I - tau * v * v'
35 : //
36 : // where tau is a real scalar, and v is a real vector with v(1:i-1) = 0
37 : // and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
38 : // and tau in TAU(i).
39 : //
40 : // The matrix P is represented in IPIV as follows: If IPIV(j) = i
41 : // then the jth column of P is the ith canonical unit vector.
42 : //
43 : // Variables in the calling sequence
44 : // ---------------------------------
45 : // M I IN The number of rows of the matrix A. M >= 0.
46 : // N I IN The number of columns of the matrix A. N >= 0.
47 : // A D IN The given M by N matrix
48 : // OUT The elements on and above the diagonal contain
49 : // the min(m,n) by n upper trapezoidal matrix R
50 : // (R is upper triangular if m >= n); the elements
51 : // belowthe diagonal, together with the array TAU,
52 : // represent the orthogonal matrix Q as a product
53 : // of min(m,n) elementary reflectors
54 : // LDA I IN The leading dimension of A. LDA >= max(1,M).
55 : // IPIV I OUT If IPIV(i) = k, then the i-th column of A*P was the
56 : // k-th column of A.
57 : // TAU I OUT Array of dimension min(M,N), The scalar factors of
58 : // the elementary reflectors
59 : // WRK D WK Work array of dimension N
60 : // SAFMIN D IN Safe minimum such that 1.0/SAFMIN does not
61 : // overflow
62 : // IER I OUT Error indicator
63 : // IER = 0 successful exit
64 : // IER = 1 data-error
65 : // IER = 2 input data error in HOUSL
66 : // IER < 0 zero pivot encountered
67 49 : int qrf(int m, int n, double *a, int lda, int *ipiv, double *tau, double *wrk,
68 : double safmin, int *ier)
69 : {
70 49 : const double fact = 0.05, one = 1.0, zer = 0.0;
71 49 : int imax = 0, ip1 = 0, mn = 0;
72 49 : double cnrm = 0.0, cnrmj = 0.0, taui = 0.0, tmp1 = 0.0, tmp2 = 0.0,
73 49 : tmp3 = 0.0;
74 :
75 : // No error has occured
76 49 : *ier = 0;
77 :
78 : // Validate input
79 49 : if ((m < 0) || (n < 0) || (lda < m))
80 : {
81 3 : *ier = 1;
82 3 : return 1;
83 : }
84 :
85 : // Quick return
86 46 : if ((m == 0) || (n == 0))
87 : {
88 2 : *ier = 0;
89 2 : return 0;
90 : }
91 :
92 : // Initialize column norms and pivot array
93 162 : for (int j = 1; j <= n; j++)
94 : {
95 118 : cnrm = dnrm2(m, a + map(1, j), 1);
96 118 : tau(j) = cnrm;
97 118 : wrk(j) = cnrm;
98 118 : ipiv(j) = j;
99 : }
100 :
101 44 : mn = n;
102 44 : if (m < n) mn = m;
103 :
104 : // Main loop for the factorization
105 156 : for (int i = 1; i <= mn; i++)
106 : {
107 116 : ip1 = i + 1;
108 : // determine pivot column
109 116 : cnrm = tau(i);
110 116 : imax = i;
111 116 : if (i < n)
112 : {
113 192 : for (int j = ip1; j <= n; j++)
114 : {
115 119 : if (tau(j) > cnrm)
116 : {
117 52 : imax = j;
118 52 : cnrm = tau(j);
119 : }
120 : }
121 73 : if (cnrm == zer)
122 : {
123 1 : printf("Message from QRF: Zero pivot encountered\n");
124 1 : *ier = -i;
125 1 : return *ier;
126 : }
127 : // Swap the columns if necessary
128 72 : if (imax != i)
129 : {
130 188 : for (int j = 1; j <= m; j++)
131 : {
132 153 : swap(double, a(j, imax), a(j, i));
133 : }
134 35 : swap(int, ipiv(imax), ipiv(i));
135 35 : tau(imax) = tau(i);
136 35 : wrk(imax) = wrk(i);
137 : }
138 : }
139 115 : if (cnrm == zer)
140 : {
141 3 : printf("Message from QRF: Zero pivot encountered\n");
142 3 : *ier = -i;
143 3 : return *ier;
144 : }
145 : else
146 : {
147 : // Generate elementary reflector H(i)
148 112 : taui = zer;
149 112 : if (i < m)
150 95 : housg(m - i + 1, a + map(i, i), a + map(ip1, i), 1, &taui, safmin);
151 112 : if (i < n)
152 : {
153 : // Apply h(i) to a(i:m,i+1:n) from the left
154 72 : housl(m - i + 1, n - i, a + map(i, i), 1, taui, a + map(i, ip1), lda,
155 : ier);
156 72 : if (*ier != 0)
157 : {
158 0 : printf("Message from QRF: Input-dimension error in housl\n");
159 0 : *ier = 2;
160 0 : return *ier;
161 : }
162 : // Update column norms
163 189 : for (int j = ip1; j <= n; j++)
164 : {
165 117 : cnrmj = tau(j);
166 117 : if (cnrmj != zer)
167 : {
168 117 : tmp1 = fabs(a(i, j)) / cnrmj;
169 117 : tmp2 = one - tmp1 * tmp1;
170 117 : if (tmp2 <= zer)
171 : {
172 2 : tmp2 = zer;
173 2 : tmp3 = zer;
174 : }
175 : else
176 : {
177 115 : tmp3 = sqrt(tmp2);
178 : }
179 117 : tmp1 = cnrmj / wrk(j);
180 117 : tmp2 = one + fact * tmp2 * tmp1 * tmp1;
181 117 : if (tmp2 == one)
182 : {
183 6 : tau(j) = dnrm2(m - i, a + map(ip1, j), 1);
184 6 : wrk(j) = tau(j);
185 : }
186 : else
187 : {
188 111 : tau(j) *= tmp3;
189 : }
190 : }
191 : }
192 : }
193 112 : tau(i) = taui;
194 : }
195 : }
196 40 : *ier = 0;
197 40 : return *ier;
198 : }
199 :
200 : // For an M x N matrix A with M >= N and rank A = N, and a given
201 : // M-vector Y, compute the least squares solution
202 : //
203 : // min { || A*X - Y || ; X in R^N }
204 : //
205 : // by the algorithm
206 : //
207 : // 1. Y := Q^T Y
208 : // 2. IF( M >= N ) THEN Solve R*Z := (I, 0)Y, Set X := Z
209 : // ELSE Solve R*Z := Y, Set X = (Z,0)
210 : //
211 : // under the assumption that the arrays A and TAU contain the
212 : // QR-factorization of an M x N matrix A.
213 : //
214 : // The routine returns Q^T Y in the array Y and the least squares
215 : // solution in the array X. The array X may be identified with Y
216 : // if Q^T*Y is not needed.
217 : //
218 : // Variables in the calling sequence
219 : // ---------------------------------
220 : // M I IN Number of rows of the matrix A, M >= N
221 : // N I IN Number of columns of the matrix A
222 : // A D IN Array of dimension LDA x n, the factored matrix
223 : // LDA I IN Leading dimension of A, LDA >= M
224 : // TAU D IN Array of dimension N containing the scalar
225 : // factor of the elementary reflectors, as
226 : // returned by QRF
227 : // Y D IN Array of dimension M, the given vector Y
228 : // OUT The vector Q^T Y
229 : // X D OUT Array of dimension N, the computed solution
230 : // IER I OUT Error indicator
231 : // IER = 0 no error
232 : // IER = 1 input data error
233 : // IER < 0 zero pivot encountered
234 : // IER = -J signifies that the
235 : // J-th diagonal element of the
236 : // triangular matrix R is zero.
237 : //
238 47 : int qrs(int m, int n, double *a, int lda, double *tau, double *y, double *x,
239 : int *ier)
240 : {
241 47 : const double zer = 0.0;
242 47 : int jm1 = 0, jp1 = 0;
243 47 : double sum = 0, t = 0;
244 :
245 : // Validate input: lda >= m >= 0 holds from these conditions
246 47 : if ((m < 0) || (n < 0) || (m < n) || (lda < m))
247 : {
248 4 : *ier = 1;
249 4 : return 1;
250 : }
251 :
252 43 : *ier = 0;
253 :
254 : // Quick return: 0 = m >= n >= 0 implies m=n=0 so that lda >= 0.
255 : // If lda=0, then there is a quick return and matrix a is not accessed.
256 43 : if ((n == 0) || ((m == 0) && (n == 0)))
257 : {
258 2 : *ier = 0;
259 2 : return 0;
260 : }
261 :
262 : // The case m = 1, under the assumption m >= n, implies n=1
263 41 : if (m == 1)
264 : {
265 2 : if (a(1, 1) == zer)
266 : {
267 1 : printf("Message from QRS: Zero pivot encountered\n");
268 1 : *ier = -1;
269 1 : return -1;
270 : }
271 : else
272 : {
273 1 : x(1) = y(1) / a(1, 1);
274 1 : *ier = 0;
275 1 : return 0;
276 : }
277 : }
278 :
279 : // Compute transpose(Q)*y
280 : //
281 : // The matrix Q is represented as a product of elementary reflectors
282 : //
283 : // Q = H(1) H(2) . . . H(k), where k = min(m,n).
284 : //
285 : // Each H(i) has the form
286 : //
287 : // H(i) = I - tau * v * v'
288 : //
289 : // v(1) = 1. transpose(Q) = H(j) H(j-1) . . . H(1).
290 147 : for (int j = 1; j <= n; j++)
291 : {
292 108 : jp1 = j + 1;
293 :
294 : // Multiply transpose(v)*y
295 108 : sum = y(j); // The first element of the elementary vector is 1
296 108 : if (j < m)
297 : {
298 437 : for (int i = jp1; i <= m; i++)
299 : {
300 346 : sum += a(i, j) * y(i); // v stored in lower triangle of A
301 : }
302 : }
303 :
304 : // y overwritten by y - tau(j)*(transpose(v)*y)*v
305 108 : if (sum != zer) // v is not zero
306 : {
307 99 : t = -tau(j) * sum;
308 99 : y(j) += t; // The first element of v is 1
309 99 : if (j < m)
310 : {
311 419 : for (int i = jp1; i <= m; i++)
312 : {
313 334 : y(i) += t * a(i, j);
314 : }
315 : }
316 : }
317 : }
318 :
319 : // Set x := the first n components of y
320 147 : for (int j = 1; j <= n; j++)
321 : {
322 108 : x(j) = y(j);
323 : }
324 : // R is stored in the upper triangle of A including diagonal.
325 : //
326 : // Solve R*Z := (I,0)Y, and set X := Z
327 : // or solve R*Z := Y, and set X := (Z,0)
328 144 : for (int j = n; j >= 1; j--)
329 : {
330 106 : if (a(j, j) == zer)
331 : {
332 1 : printf("Message from QRS: Zero pivot encountered\n");
333 1 : *ier = -j;
334 1 : return -j;
335 : }
336 : // Back substitution
337 105 : x(j) = x(j) / a(j, j);
338 105 : if (j > 1)
339 : {
340 67 : jm1 = j - 1;
341 67 : t = -x(j);
342 176 : for (int i = 1; i <= jm1; i++)
343 : {
344 109 : x(i) += t * a(i, j);
345 : }
346 : }
347 : }
348 38 : *ier = 0;
349 38 : return *ier;
350 : }
351 :
352 : // Generate an M by NQ real matrix Q with NQ = L - K + 1 orthonormal
353 : // columns formed as the product of N Householder reflectors
354 : // H(1),...,H(N) of order M >= N as returned in the columns of the
355 : // array A by QRF.
356 : //
357 : // Variables in the calling sequence:
358 : // ----------------------------------
359 : //
360 : // M I IN The number of rows of the matrix A. M >= 0.
361 : // N I IN The number of columns of the matrix A. M >= N >= 0.
362 : // K I IN The lower column index, M >= K >= 1
363 : // K is set to 1 if K <= 0
364 : // L I IN The higher column index, M >= L >= K >= 1
365 : // L is set to M if L >= M
366 : // A D IN Array of dimension (LDA,N)
367 : // For IA = 1 the i-th column of A is assumed to contain
368 : // the vector defining the elementary reflector
369 : // H(i), for i = 1,2,...,N, as returned by QRF
370 : // For IA = 2 the i-th row of A is assumed to contain
371 : // the vector defining the elementary reflector
372 : // H(i), for i = 1,2,...,M, as returned by QRF
373 : // LDA I IN The leading dimension of the array A, LDA >= max(1,M).
374 : // TAU D IN Array of dimension (K) where TAU(i) contains the
375 : // scalar factor of the elementary reflector H(i), as
376 : // returned by QRF.
377 : // Q D OUT The orthogonal matrix.
378 : // LDQ D IN The leading dimension of the array Q, LDQ >= max(1,M)
379 : // IER I OUT Error indicator
380 : // IER = 0 successful exit
381 : // IER = 1 input-data error
382 : // IER = 2 input-data error in HOUSL
383 : //
384 : // NOTE: The arrays A and Q cannot be identified
385 9 : int qorg(int m, int n, int k, int l, double *a, int lda, double *tau, double *q,
386 : int ldq, int *ier)
387 : {
388 9 : const double zer = 0.0, one = 1.0;
389 9 : int j2 = 0, km1 = 0, nl = 0, nq = 0;
390 :
391 : // Validate input: lda >= m >= 0 holds from these conditions
392 9 : if ((m < 0) || (n < 0) || (m < n) || (k > l) || (lda < m) || (ldq < 1) ||
393 : (ldq < m))
394 : {
395 6 : *ier = 1;
396 6 : return 1;
397 : }
398 :
399 3 : *ier = 0;
400 :
401 : // Quick return: 0 = m >= n >= 0 implies m=n=0 so that lda >= 0.
402 : // If lda=0, then there is a quick return and matrix a is not accessed.
403 : // k=0 handled because k < = 0 means k set to 1 below.
404 3 : if ((n == 0) || ((m == 0) && (n == 0)))
405 : {
406 2 : *ier = 0;
407 2 : return 0;
408 : }
409 :
410 : // Set defaults
411 1 : if (k < 0) k = 1;
412 1 : if (l > m) l = m;
413 :
414 : // Initialize Q
415 1 : km1 = k - 1;
416 1 : nq = l - km1;
417 4 : for (int j1 = 1; j1 <= nq; j1++)
418 : {
419 12 : for (int i = 1; i <= m; i++)
420 : {
421 9 : q(i, j1) = zer;
422 : }
423 3 : q(j1 + km1, j1) = one;
424 : }
425 :
426 : // Apply H(1)*...*H(NL)*Q
427 1 : nl = l;
428 1 : if (l > n) nl = l;
429 1 : j2 = nq;
430 4 : for (int j1 = nl; j1 >= 1; j1--)
431 : {
432 : // Apply H(J1) to Q(J1:M,J1:NQ) from the left
433 3 : housl(m - j1 + 1, nq - j2 + 1, a + map(j1, j1), 1, tau(j1), q + map(j1, j2),
434 : ldq, ier);
435 3 : if (*ier != 0)
436 : {
437 0 : printf("Message from QORG: Input-dimension error in housl\n");
438 0 : *ier = 2;
439 0 : return *ier;
440 : }
441 3 : j2 = j2 - 1;
442 3 : if (j2 < 1) j2 = 1;
443 : }
444 1 : *ier = 0;
445 1 : return *ier;
446 : }
447 :
448 : // Generates an n-dimensional Householder reflector
449 : //
450 : // H = I - tau*( 1 ) * ( 1 v' ), H' * H = I, tau scalar
451 : // ( v )
452 : //
453 : // such that
454 : //
455 : // H * ( alpha ) = ( beta ), alpha, beta scalars
456 : // ( x ) ( 0 ) x (n-1)-dimensional vector
457 : //
458 : // Because of H'* H = I it follows that
459 : //
460 : // alpha^2 + x^T x = beta^2 ==> beta = sqrt(alpha^2 + x^T x)
461 : //
462 : // H'( beta ) = ( alpha ) ==> tau = (beta - alpha)/beta
463 : // ( 0 ) ( x ) v = xscal*x,
464 : // xscal = 1/(beta - alpha)
465 : //
466 : // If x = 0, then tau = 0 and H = I, otherwise 1 <= tau <= 2.
467 : //
468 : // This is an edited version of the LAPACK routine DLARFG
469 : //
470 : // Variables in the calling sequence:
471 : // ----------------------------------
472 : // N I IN Dimension of H
473 : // ALPHA D IN The scalar alpha
474 : // OUT The scalar beta
475 : // X D IN The given vector x of dimension n - 1
476 : // OUT The vector v
477 : // INCX I IN The increment between elements of X, INCX .NE. 0
478 : // TAU D OUT The scalar tau
479 : // SAFMIN D IN Safe minimum such that 1.0/SAFMIN does not
480 : // overflow
481 106 : int housg(int n, double *alpha, double *x, int incx, double *tau, double safmin)
482 : {
483 106 : const double one = 1.0, zer = 0.0;
484 106 : int ix = 0, knt = 0, kx = 0, nm1 = 0;
485 106 : double a1 = 0.0, a2 = 0.0, beta = 0.0, rsafmn = 0.0, tmp = 0.0, xnorm = 0.0,
486 106 : xscal = 0.0;
487 :
488 : // H is defined as the identity whenever incx==0; implies don't iterate
489 : // through x
490 106 : if (incx == 0)
491 : {
492 1 : *tau = zer;
493 1 : return -2;
494 : }
495 :
496 : // H is defined as the identity whenever n<1
497 105 : if (n < 1)
498 : {
499 1 : *tau = zer;
500 1 : return -1;
501 : }
502 :
503 : // H is defined as the identity whenever n=1
504 104 : if (n == 1)
505 : {
506 1 : *tau = zer;
507 1 : return 0;
508 : }
509 :
510 : // Dimension of x
511 103 : nm1 = n - 1;
512 : // Norm of x
513 103 : xnorm = dnrm2(nm1, x, incx);
514 :
515 : // H is the identity whenever xnorm=0 with alpha=beta and tau=0
516 103 : if (xnorm == zer)
517 : {
518 5 : *tau = zer;
519 5 : return 0;
520 : }
521 :
522 : // General case
523 98 : kx = 1;
524 98 : if (incx < 0) kx = 1 - (n - 1) * incx;
525 :
526 : // Compute alpha, beta, tau, and xscal
527 98 : a1 = xnorm;
528 98 : a2 = fabs(*alpha);
529 98 : if (a1 < a2)
530 : {
531 45 : a1 = a2;
532 45 : a2 = xnorm;
533 : }
534 98 : if (a2 == zer)
535 : {
536 2 : beta = a1;
537 : }
538 : else
539 : {
540 96 : tmp = a1 / a2;
541 96 : beta = a2 * sqrt(one + tmp * tmp);
542 : }
543 98 : if (*alpha > zer) beta = -beta;
544 :
545 : // Test for loss of accuracy
546 98 : if (fabs(beta) >= safmin)
547 : {
548 98 : *tau = (beta - *alpha) / beta;
549 98 : xscal = one / (*alpha - beta);
550 98 : *alpha = beta;
551 : }
552 : else // xnorm, beta may be inaccurate; scale x and recompute
553 : {
554 0 : rsafmn = one / safmin;
555 0 : knt = 0;
556 : do
557 : {
558 0 : knt = knt + 1;
559 0 : if (incx == 1)
560 : {
561 0 : for (int i = 1; i <= nm1; i++)
562 : {
563 0 : x(i) = rsafmn * x(i);
564 : }
565 : }
566 : else
567 : {
568 0 : ix = kx;
569 0 : for (int i = 1; i <= nm1; i++)
570 : {
571 0 : x(ix) = rsafmn * x(ix);
572 0 : ix = ix + incx;
573 : }
574 : }
575 0 : beta = beta * rsafmn;
576 0 : *alpha = (*alpha) * rsafmn;
577 0 : } while (fabs(beta) < safmin);
578 :
579 : // new beta satisfies safmin <= beta <= 1.0
580 0 : xnorm = dnrm2(nm1, x, incx);
581 0 : a1 = xnorm;
582 0 : a2 = fabs(*alpha);
583 0 : if (a1 < a2)
584 : {
585 0 : a1 = a2;
586 0 : a2 = xnorm;
587 : }
588 0 : if (a2 == zer)
589 : {
590 0 : beta = a1;
591 : }
592 : else
593 : {
594 0 : tmp = a1 / a2;
595 0 : beta = a2 * sqrt(one + tmp * tmp);
596 : }
597 0 : if (*alpha > zer) beta = -beta;
598 0 : *tau = (beta - *alpha) / beta;
599 0 : xscal = one / (*alpha - beta);
600 0 : *alpha = beta;
601 0 : for (int j = 1; j <= knt; j++)
602 : {
603 0 : *alpha = (*alpha) * safmin;
604 : }
605 : }
606 98 : if (incx == 1)
607 : {
608 490 : for (int i = 1; i <= nm1; i++)
609 : {
610 392 : x(i) = xscal * x(i);
611 : }
612 : }
613 : else
614 : {
615 0 : ix = kx;
616 0 : for (int i = 1; i <= nm1; i++)
617 : {
618 0 : x(ix) = xscal * x(ix);
619 0 : ix = ix + incx;
620 : }
621 : }
622 98 : return 0;
623 : }
624 :
625 : // Multiplies a given M x N matrix A from the (L)eft by an
626 : // M-dimensional Householder reflector
627 : //
628 : // H = I - tau*x*x', H' * H = I, x = ( 1 )
629 : // ( v )
630 : //
631 : // that is, overwrite A by the product H * A
632 : //
633 : // Variables in the calling sequence:
634 : // ----------------------------------
635 : // M I IN The number of rows of A
636 : // N I IN The number of columns of A
637 : // X D IN The given vector x. x(1) = 1.0 is enforced
638 : // INCX I IN The increment between elements of X, INCX .NE. 0
639 : // TAU D IN The scalar tau
640 : // A D IN The given matrix of dimension M x N
641 : // OUT The computed matrix product H * A
642 : // LDA I IN The leading dimension of the array A, LDA >= M
643 : // IER I OUT Error indicator
644 : // IER = 0 no error
645 : // IER = 1 input-data error
646 85 : int housl(int m, int n, double *x, int incx, double tau, double *a, int lda,
647 : int *ier)
648 : {
649 85 : const double zer = 0.0;
650 85 : int ix = 0, kx = 0;
651 85 : double sum = 0.0, tmp = 0.0;
652 :
653 85 : *ier = 0;
654 :
655 : // Valid data
656 85 : if ((m < 0) || (n < 0) || (incx == 0) || (lda < m))
657 : {
658 4 : *ier = 1;
659 4 : return 1;
660 : }
661 :
662 : // Quick return
663 81 : if ((m == 0) || (n == 0) || (tau == 0.0))
664 : {
665 7 : *ier = 0;
666 7 : return 0;
667 : }
668 :
669 74 : kx = 1;
670 74 : if (incx < 0) kx = 1 - (m - 1) * incx;
671 :
672 : // Compute w = A' * x and A := A - tau * x * w'
673 :
674 200 : for (int j = 1; j <= n; j++)
675 : {
676 126 : sum = a(1, j);
677 126 : if (m > 1)
678 : {
679 126 : ix = kx;
680 506 : for (int i = 2; i <= m; i++)
681 : {
682 380 : ix += incx;
683 380 : sum += a(i, j) * x(ix);
684 : }
685 : }
686 126 : if (sum != zer)
687 : {
688 120 : tmp = -tau * sum;
689 120 : a(1, j) += tmp;
690 120 : if (m > 1)
691 : {
692 120 : ix = kx;
693 484 : for (int i = 2; i <= m; i++)
694 : {
695 364 : ix += incx;
696 364 : a(i, j) += x(ix) * tmp;
697 : }
698 : }
699 : }
700 : }
701 74 : return *ier;
702 : }
|
