1 package ubic.basecode.math.linalg;
2
3 /*
4 QR_j.java copyright claim:
5
6 This software is based on public domain LINPACK routines.
7 It was translated from FORTRAN to Java by a US government employee
8 on official time. Thus this software is also in the public domain.
9
10
11 The translator's mail address is:
12
13 Steve Verrill
14 USDA Forest Products Laboratory
15 1 Gifford Pinchot Drive
16 Madison, Wisconsin
17 53705
18
19
20 The translator's e-mail address is:
21
22 steve@swst.org
23
24
25 ***********************************************************************
26
27 DISCLAIMER OF WARRANTIES:
28
29 THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND.
30 THE TRANSLATOR DOES NOT WARRANT, GUARANTEE OR MAKE ANY REPRESENTATIONS
31 REGARDING THE SOFTWARE OR DOCUMENTATION IN TERMS OF THEIR CORRECTNESS,
32 RELIABILITY, CURRENTNESS, OR OTHERWISE. THE ENTIRE RISK AS TO
33 THE RESULTS AND PERFORMANCE OF THE SOFTWARE IS ASSUMED BY YOU.
34 IN NO CASE WILL ANY PARTY INVOLVED WITH THE CREATION OR DISTRIBUTION
35 OF THE SOFTWARE BE LIABLE FOR ANY DAMAGE THAT MAY RESULT FROM THE USE
36 OF THIS SOFTWARE.
37
38 Sorry about that.
39
40 ***********************************************************************
41
42
43 History:
44
45 Date Translator Changes
46
47 2/25/97 Steve Verrill Translated
48 6/5/97 Java/C style indexing
49
50 */
51
52 /**
53 *
54 *
55 *
56 * This class contains the LINPACK DQRDC (QR decomposition)
57 * and DQRSL (QR solve) routines.
58 *
59 *
60 *
61 *
62 * IMPORTANT: The "_j" suffixes indicate that these routines use
63 * Java/C style indexing. For example, you will see
64 *
65 * for (i = 0; i < n; i++)
66 *
67 * rather than
68 *
69 * for (i = 1; i <= n; i++)
70 *
71 * To use the "_j" routines you will have
72 * to fill elements 0 through n - 1 rather than elements 1 through n.
73 * Versions of these programs that use FORTYRAN style indexing are
74 * also available. They end with the suffix "_f77".
75 *
76 *
77 *
78 *
79 * This class was translated by a statistician from FORTRAN versions of
80 * the LINPACK routines. It is NOT an official translation. When
81 * public domain Java numerical analysis routines become available
82 * from the people who produce LAPACK, then THE CODE PRODUCED
83 * BY THE NUMERICAL ANALYSTS SHOULD BE USED.
84 *
85 *
86 *
87 *
88 * Meanwhile, if you have suggestions for improving this
89 * code, please contact Steve Verrill at sverrill@fs.fed.us.
90 *
91 * @author Steve Verrill
92 * @version .5 --- June 5, 1997
93 *
94 */
95 public class Dqrsl extends Object {
96
97 /**
98 *
99 *
100 * This method decomposes an n by p matrix X into a product, QR, of
101 * an orthogonal n by n matrix Q and an upper triangular n by p matrix R.
102 * For details, see the comments in the code.
103 * This method is a translation from FORTRAN to Java of the LINPACK subroutine
104 * DQRDC. In the LINPACK listing DQRDC is attributed to G.W. Stewart
105 * with a date of 8/14/78.
106 *
107 * Translated by Steve Verrill, February 25, 1997.
108 *
109 * @param X The matrix to be decomposed
110 * @param n The number of rows of the matrix X
111 * @param p The number of columns of the matrix X
112 * @param qraux This vector "contains further information required to
113 * recover the orthogonal part of the decomposition."
114 * @param jpvt This output vector contains pivoting information.
115 * @param work This vector is used as temporary space
116 * @param job This value indicates whether column pivoting should be performed
117 *
118 */
119 public static void dqrdc_j( double x[][], int n, int p, double qraux[], int jpvt[], int job ) {
120
121 int j, jj, jp, l, lp1, lup, maxj, pl, pu;
122 int jm1, plm1, pum1, lm1, maxjm1;
123
124 double maxnrm, tt;
125 double nrmxl, t;
126 double fac;
127 double work[] = new double[p];
128
129 boolean initial, fin;
130
131 /*
132 *
133 * Here is a copy of the LINPACK documentation (from the SLATEC version):
134 *
135 * C***BEGIN PROLOGUE DQRDC
136 * C***DATE WRITTEN 780814 (YYMMDD)
137 * C***REVISION DATE 861211 (YYMMDD)
138 * C***CATEGORY NO. D5
139 * C***KEYWORDS LIBRARY=SLATEC(LINPACK),
140 * C TYPE=DOUBLE PRECISION(SQRDC-S DQRDC-D CQRDC-C),
141 * C LINEAR ALGEBRA,MATRIX,ORTHOGONAL TRIANGULAR,
142 * C QR DECOMPOSITION
143 * C***AUTHOR STEWART, G. W., (U. OF MARYLAND)
144 * C***PURPOSE Uses Householder transformations to compute the QR factori-
145 * C zation of N by P matrix X. Column pivoting is optional.
146 * C***DESCRIPTION
147 * C
148 * C DQRDC uses Householder transformations to compute the QR
149 * C factorization of an N by P matrix X. Column pivoting
150 * C based on the 2-norms of the reduced columns may be
151 * C performed at the user's option.
152 * C
153 * C On Entry
154 * C
155 * C X DOUBLE PRECISION(LDX,P), where LDX .GE. N.
156 * C X contains the matrix whose decomposition is to be
157 * C computed.
158 * C
159 * C LDX INTEGER.
160 * C LDX is the leading dimension of the array X.
161 * C
162 * C N INTEGER.
163 * C N is the number of rows of the matrix X.
164 * C
165 * C P INTEGER.
166 * C P is the number of columns of the matrix X.
167 * C
168 * C JPVT INTEGER(P).
169 * C JPVT contains integers that control the selection
170 * C of the pivot columns. The K-th column X(K) of X
171 * C is placed in one of three classes according to the
172 * C value of JPVT(K).
173 * C
174 * C If JPVT(K) .GT. 0, then X(K) is an initial
175 * C column.
176 * C
177 * C If JPVT(K) .EQ. 0, then X(K) is a free column.
178 * C
179 * C If JPVT(K) .LT. 0, then X(K) is a final column.
180 * C
181 * C Before the decomposition is computed, initial columns
182 * C are moved to the beginning of the array X and final
183 * C columns to the end. Both initial and final columns
184 * C are frozen in place during the computation and only
185 * C free columns are moved. At the K-th stage of the
186 * C reduction, if X(K) is occupied by a free column
187 * C it is interchanged with the free column of largest
188 * C reduced norm. JPVT is not referenced if
189 * C JOB .EQ. 0.
190 * C
191 * C WORK DOUBLE PRECISION(P).
192 * C WORK is a work array. WORK is not referenced if
193 * C JOB .EQ. 0.
194 * C
195 * C JOB INTEGER.
196 * C JOB is an integer that initiates column pivoting.
197 * C If JOB .EQ. 0, no pivoting is done.
198 * C If JOB .NE. 0, pivoting is done.
199 * C
200 * C On Return
201 * C
202 * C X X contains in its upper triangle the upper
203 * C triangular matrix R of the QR factorization.
204 * C Below its diagonal X contains information from
205 * C which the orthogonal part of the decomposition
206 * C can be recovered. Note that if pivoting has
207 * C been requested, the decomposition is not that
208 * C of the original matrix X but that of X
209 * C with its columns permuted as described by JPVT.
210 * C
211 * C QRAUX DOUBLE PRECISION(P).
212 * C QRAUX contains further information required to recover
213 * C the orthogonal part of the decomposition.
214 * C
215 * C JPVT JPVT(K) contains the index of the column of the
216 * C original matrix that has been interchanged into
217 * C the K-th column, if pivoting was requested.
218 * C
219 * C LINPACK. This version dated 08/14/78 .
220 * C G. W. Stewart, University of Maryland, Argonne National Lab.
221 * C
222 * C DQRDC uses the following functions and subprograms.
223 * C
224 * C BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2
225 * C Fortran DABS,DMAX1,MIN0,DSQRT
226 * C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,
227 * C *LINPACK USERS GUIDE*, SIAM, 1979.
228 * C***ROUTINES CALLED DAXPY,DDOT,DNRM2,DSCAL,DSWAP
229 * C***END PROLOGUE DQRDC
230 *
231 */
232
233 pl = 1;
234 plm1 = pl - 1;
235 pu = 0;
236
237 if ( job != 0 ) {
238
239 // Pivoting has been requested. Rearrange the columns according to pvt.
240
241 for ( j = 1; j <= p; j++ ) {
242 jm1 = j - 1;
243 initial = ( jpvt[jm1] > 0 );
244 fin = ( jpvt[jm1] < 0 );
245
246 jpvt[jm1] = j;
247 if ( fin ) jpvt[jm1] = -j;
248
249 if ( initial ) {
250 if ( j != pl ) Blas.colswap_j( n, x, plm1, jm1 );
251 jpvt[jm1] = jpvt[plm1];
252 jpvt[plm1] = j;
253 plm1 = pl;
254 pl++;
255 }
256 }
257
258 pu = p;
259 pum1 = p - 1;
260
261 for ( jj = 1; jj <= p; jj++ ) {
262
263 j = p - jj;
264
265 if ( jpvt[j] < 0 ) {
266
267 jpvt[j] = -jpvt[j];
268
269 if ( j != pum1 ) {
270 Blas.colswap_j( n, x, pum1, j );
271 jp = jpvt[pum1];
272 jpvt[pum1] = jpvt[j];
273 jpvt[j] = jp;
274 }
275
276 pu--;
277 pum1--;
278
279 }
280 }
281 }
282
283 // Compute the norms of the free columns.
284 for ( j = pl - 1; j < pu; j++ ) {
285 qraux[j] = Blas.colnrm2_j( n, x, 0, j );
286 work[j] = qraux[j];
287 }
288
289 // Perform the Householder reduction of X.
290
291 lup = Math.min( n, p );
292
293 for ( l = 1; l <= lup; l++ ) {
294
295 lm1 = l - 1;
296
297 if ( l >= pl && l < pu ) {
298
299 // Locate the column of greatest norm and bring it into
300 // the pivot position.
301
302 maxnrm = 0.0;
303 maxj = l;
304
305 for ( j = l; j <= pu; j++ ) {
306
307 jm1 = j - 1;
308
309 if ( qraux[jm1] > maxnrm ) {
310 maxnrm = qraux[jm1];
311 maxj = j;
312 }
313 }
314
315 if ( maxj != l ) {
316 maxjm1 = maxj - 1;
317 Blas.colswap_j( n, x, lm1, maxjm1 );
318 qraux[maxjm1] = qraux[lm1];
319 work[maxjm1] = work[lm1];
320 jp = jpvt[maxjm1];
321 jpvt[maxjm1] = jpvt[lm1];
322 jpvt[lm1] = jp;
323
324 }
325 }
326
327 qraux[lm1] = 0.0;
328
329 if ( l != n ) {
330
331 // Compute the Householder transformation for column l.
332
333 nrmxl = Blas.colnrm2_j( n - l + 1, x, lm1, lm1 );
334
335 if ( nrmxl != 0.0 ) {
336
337 if ( x[lm1][lm1] != 0.0 ) nrmxl = Blas.sign_j( nrmxl, x[lm1][lm1] );
338
339 Blas.colscal_j( n - l + 1, 1.0 / nrmxl, x, lm1, lm1 );
340
341 x[lm1][lm1]++;
342
343 // Apply the transformation to the remaining columns,
344 // updating the norms.
345
346 lp1 = l + 1;
347
348 for ( j = lp1; j <= p; j++ ) {
349
350 jm1 = j - 1;
351
352 t = -Blas.coldot_j( n - l + 1, x, lm1, lm1, jm1 ) / x[lm1][lm1];
353
354 Blas.colaxpy_j( n - l + 1, t, x, lm1, lm1, jm1 );
355
356 if ( j >= pl && j <= pu ) {
357
358 if ( qraux[jm1] != 0.0 ) {
359
360 fac = Math.abs( x[lm1][jm1] ) / qraux[jm1];
361 tt = 1.0 - fac * fac;
362 tt = Math.max( tt, 0.0 );
363 t = tt;
364 fac = qraux[jm1] / work[jm1];
365 tt = 1.0 + .05 * tt * fac * fac;
366
367 if ( tt != 1.0 ) {
368
369 qraux[jm1] *= Math.sqrt( t );
370
371 } else {
372
373 qraux[jm1] = Blas.colnrm2_j( n - l, x, l, jm1 );
374 work[jm1] = qraux[jm1];
375
376 }
377 }
378 }
379 }
380
381 // Save the transformation
382
383 qraux[lm1] = x[lm1][lm1];
384 x[lm1][lm1] = -nrmxl;
385
386 }
387 }
388 }
389
390 return;
391
392 }
393
394 /**
395 * This method "applies the output of DQRDC to compute coordinate
396 * transformations, projections, and least squares solutions."
397 * For details, see the comments in the code.
398 * This method is a translation from FORTRAN to Java of the LINPACK subroutine
399 * DQRSL. In the LINPACK listing DQRSL is attributed to G.W. Stewart
400 * with a date of 8/14/78.
401 *
402 * Translated by Steve Verrill, February 27, 1997.
403 *
404 * @param X This n by p matrix contains most of the output from DQRDC
405 * @param n The number of rows of X
406 * @param k k <= min(n,p) where p is the number of columns of X
407 * @param qraux This vector "contains further information required to
408 * recover the orthogonal part of the decomposition"
409 * @param y This n by 1 vector will be manipulated by DQRSL
410 * @param qy On output, this vector contains Qy if it has been requested
411 * @param qty On output, this vector contains transpose(Q)y if it has been requested
412 * @param b Parameter estimates
413 * @param rsd Residuals
414 * @param xb Predicted values
415 * @param job Specifies what is to be computed (see the code for details)
416 *
417 */
418 public static int dqrsl_j( double x[][], int n, int k, double qraux[],
419 double y[], double qy[], double qty[],
420 double b[], double rsd[], double xb[],
421 int job ) {
422
423 /*
424 *
425 * Here is a copy of the LINPACK documentation (from the SLATEC version):
426 *
427 * c***begin prologue dqrsl
428 * c***date written 780814 (yymmdd)
429 * c***revision date 861211 (yymmdd)
430 * c***category no. d9,d2a1
431 * c***keywords library=slatec(linpack),
432 * c type=double precision(sqrsl-s dqrsl-d cqrsl-c),
433 * c linear algebra,matrix,orthogonal triangular,solve
434 * c***author stewart, g. w., (u. of maryland)
435 * c***purpose applies the output of dqrdc to compute coordinate
436 * c transformations, projections, and least squares solutions.
437 * c***description
438 * c
439 * c dqrsl applies the output of dqrdc to compute coordinate
440 * c transformations, projections, and least squares solutions.
441 * c for k .le. min(n,p), let xk be the matrix
442 * c
443 * c xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k)))
444 * c
445 * c formed from columnns jpvt(1), ... ,jpvt(k) of the original
446 * c n x p matrix x that was input to dqrdc (if no pivoting was
447 * c done, xk consists of the first k columns of x in their
448 * c original order). dqrdc produces a factored orthogonal matrix q
449 * c and an upper triangular matrix r such that
450 * c
451 * c xk = q * (r)
452 * c (0)
453 * c
454 * c this information is contained in coded form in the arrays
455 * c x and qraux.
456 * c
457 * c on entry
458 * c
459 * c x double precision(ldx,p).
460 * c x contains the output of dqrdc.
461 * c
462 * c ldx integer.
463 * c ldx is the leading dimension of the array x.
464 * c
465 * c n integer.
466 * c n is the number of rows of the matrix xk. it must
467 * c have the same value as n in dqrdc.
468 * c
469 * c k integer.
470 * c k is the number of columns of the matrix xk. k
471 * c must not be greater than min(n,p), where p is the
472 * c same as in the calling sequence to dqrdc.
473 * c
474 * c qraux double precision(p).
475 * c qraux contains the auxiliary output from dqrdc.
476 * c
477 * c y double precision(n)
478 * c y contains an n-vector that is to be manipulated
479 * c by dqrsl.
480 * c
481 * c job integer.
482 * c job specifies what is to be computed. job has
483 * c the decimal expansion abcde, with the following
484 * c meaning.
485 * c
486 * c if a .ne. 0, compute qy.
487 * c if b,c,d, or e .ne. 0, compute qty.
488 * c if c .ne. 0, compute b.
489 * c if d .ne. 0, compute rsd.
490 * c if e .ne. 0, compute xb.
491 * c
492 * c note that a request to compute b, rsd, or xb
493 * c automatically triggers the computation of qty, for
494 * c which an array must be provided in the calling
495 * c sequence.
496 * c
497 * c on return
498 * c
499 * c qy double precision(n).
500 * c qy contains q*y, if its computation has been
501 * c requested.
502 * c
503 * c qty double precision(n).
504 * c qty contains trans(q)*y, if its computation has
505 * c been requested. here trans(q) is the
506 * c transpose of the matrix q.
507 * c
508 * c b double precision(k)
509 * c b contains the solution of the least squares problem
510 * c
511 * c minimize norm2(y - xk*b),
512 * c
513 * c if its computation has been requested. (note that
514 * c if pivoting was requested in dqrdc, the j-th
515 * c component of b will be associated with column jpvt(j)
516 * c of the original matrix x that was input into dqrdc.)
517 * c
518 * c rsd double precision(n).
519 * c rsd contains the least squares residual y - xk*b,
520 * c if its computation has been requested. rsd is
521 * c also the orthogonal projection of y onto the
522 * c orthogonal complement of the column space of xk.
523 * c
524 * c xb double precision(n).
525 * c xb contains the least squares approximation xk*b,
526 * c if its computation has been requested. xb is also
527 * c the orthogonal projection of y onto the column space
528 * c of x.
529 * c
530 * c info integer.
531 * c info is zero unless the computation of b has
532 * c been requested and r is exactly singular. in
533 * c this case, info is the index of the first zero
534 * c diagonal element of r and b is left unaltered.
535 *
536 *
537 * For the Java version, info
538 * is the return value of the the dqrsl_j method.
539 *
540 *
541 * c
542 * c the parameters qy, qty, b, rsd, and xb are not referenced
543 * c if their computation is not requested and in this case
544 * c can be replaced by dummy variables in the calling program.
545 * c to save storage, the user may in some cases use the same
546 * c array for different parameters in the calling sequence. a
547 * c frequently occuring example is when one wishes to compute
548 * c any of b, rsd, or xb and does not need y or qty. in this
549 * c case one may identify y, qty, and one of b, rsd, or xb, while
550 * c providing separate arrays for anything else that is to be
551 * c computed. thus the calling sequence
552 * c
553 * c call dqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info)
554 * c
555 * c will result in the computation of b and rsd, with rsd
556 * c overwriting y. more generally, each item in the following
557 * c list contains groups of permissible identifications for
558 * c a single calling sequence.
559 * c
560 * c 1. (y,qty,b) (rsd) (xb) (qy)
561 * c
562 * c 2. (y,qty,rsd) (b) (xb) (qy)
563 * c
564 * c 3. (y,qty,xb) (b) (rsd) (qy)
565 * c
566 * c 4. (y,qy) (qty,b) (rsd) (xb)
567 * c
568 * c 5. (y,qy) (qty,rsd) (b) (xb)
569 * c
570 * c 6. (y,qy) (qty,xb) (b) (rsd)
571 * c
572 * c in any group the value returned in the array allocated to
573 * c the group corresponds to the last member of the group.
574 * c
575 * c linpack. this version dated 08/14/78 .
576 * c g. w. stewart, university of maryland, argonne national lab.
577 * c
578 * c dqrsl uses the following functions and subprograms.
579 * c
580 * c blas daxpy,dcopy,ddot
581 * c fortran dabs,min0,mod
582 * c***references dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
583 * c *linpack users guide*, siam, 1979.
584 * c***routines called daxpy,dcopy,ddot
585 * c***end prologue dqrsl
586 *
587 */
588
589 /*
590 * Input, integer JOB, specifies what is to be computed. JOB has
591 * the decimal expansion ABCDE, with the following meaning:
592 * if A != 0, compute QY.
593 * if B != 0, compute QTY.
594 * if C != 0, compute QTY and B.
595 * if D != 0, compute QTY and RSD.
596 * if E != 0, compute QTY and AB.
597 */
598
599 int i, j, jj, ju, info;
600 int jm1;
601 double t, temp;
602 boolean cb, cqy, cqty, cr, cxb;
603
604 // set info flag
605
606 info = 0;
607
608 // determine what is to be computed
609
610 cqy = ( job / 10000 != 0 );
611 cqty = ( ( job % 10000 ) != 0 );
612 cb = ( ( job % 1000 ) / 100 != 0 );
613 cr = ( ( job % 100 ) / 10 != 0 );
614 cxb = ( ( job % 10 ) != 0 );
615
616 ju = Math.min( k, n - 1 );
617
618 // special action when n = 1
619
620 if ( ju == 0 ) {
621
622 if ( cqy ) qy[0] = y[0];
623 if ( cqty ) qty[0] = y[0];
624 if ( cxb ) xb[0] = y[0];
625
626 if ( cb ) {
627
628 if ( x[0][0] == 0.0 ) {
629
630 info = 1;
631
632 } else {
633
634 b[0] = y[0] / x[0][0];
635
636 }
637
638 }
639
640 if ( cr ) rsd[0] = 0.0;
641
642 return info;
643
644 }
645
646 // set up to compute qy or transpose(q)y
647
648 if ( cqy ) Blas.dcopy_j( n, y, 1, qy, 1 );
649 if ( cqty ) Blas.dcopy_j( n, y, 1, qty, 1 );
650
651 if ( cqy ) {
652
653 // compute qy
654
655 for ( jj = 1; jj <= ju; jj++ ) {
656
657 j = ju - jj;
658
659 if ( qraux[j] != 0.0 ) {
660
661 temp = x[j][j];
662 x[j][j] = qraux[j];
663 t = -Blas.colvdot_j( n - j, x, qy, j, j ) / x[j][j];
664 Blas.colvaxpy_j( n - j, t, x, qy, j, j );
665 x[j][j] = temp;
666
667 }
668
669 }
670
671 }
672
673 if ( cqty ) {
674
675 // compute transpose(q)y
676
677 for ( j = 0; j < ju; j++ ) {
678 if ( qraux[j] != 0.0 ) {
679 temp = x[j][j];
680 x[j][j] = qraux[j];
681 t = -Blas.colvdot_j( n - j, x, qty, j, j ) / x[j][j];
682 Blas.colvaxpy_j( n - j, t, x, qty, j, j );
683 x[j][j] = temp;
684 }
685 }
686 }
687
688 // set up to compute b, rsd, or xb
689
690 if ( cb ) Blas.dcopy_j( k, qty, 1, b, 1 );
691
692 if ( cxb ) Blas.dcopy_j( k, qty, 1, xb, 1 );
693
694 if ( cr && ( k < n ) ) Blas.dcopyp_j( n - k, qty, rsd, k );
695
696 if ( cxb ) {
697 for ( i = k; i < n; i++ ) {
698 xb[i] = 0.0;
699 }
700 }
701
702 if ( cr ) {
703 for ( i = 0; i < k; i++ ) {
704 rsd[i] = 0.0;
705 }
706 }
707
708 if ( cb ) {
709
710 // compute b
711
712 for ( jj = 1; jj <= k; jj++ ) {
713
714 jm1 = k - jj;
715 if ( x[jm1][jm1] == 0.0 ) {
716 info = jm1 + 1;
717 break;
718 }
719
720 b[jm1] = b[jm1] / x[jm1][jm1];
721
722 if ( jm1 != 0 ) {
723 t = -b[jm1];
724 Blas.colvaxpy_j( jm1, t, x, b, 0, jm1 );
725
726 }
727 }
728 }
729
730 if ( cr || cxb ) {
731 // compute rsd or xb as required
732 for ( jj = 1; jj <= ju; jj++ ) {
733 j = ju - jj;
734
735 if ( qraux[j] != 0.0 ) {
736 temp = x[j][j];
737 x[j][j] = qraux[j];
738
739 if ( cr ) {
740 t = -Blas.colvdot_j( n - j, x, rsd, j, j ) / x[j][j];
741 Blas.colvaxpy_j( n - j, t, x, rsd, j, j );
742 }
743
744 if ( cxb ) {
745 t = -Blas.colvdot_j( n - j, x, xb, j, j ) / x[j][j];
746 Blas.colvaxpy_j( n - j, t, x, xb, j, j );
747
748 }
749
750 x[j][j] = temp;
751
752 }
753 }
754
755 }
756
757 return info;
758
759 }
760
761 }