1 /*
2 *
3 * https://www1.fpl.fs.fed.us/linear_algebra.html
4 Blas_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/3/97 Steve Verrill Java/C style indexing
49 3/10/98 Steve Verrill isamax and colisamax added
50
51 */
52
53 package ubic.basecode.math.linalg;
54
55 /**
56 *
57 * <p>
58 * This class contains Java versions of a number of the LINPACK
59 * basic linear algebra subroutines (blas):
60 * <ol>
61 * <li>isamax_j
62 * <li>daxpy_j
63 * <li>ddot_j
64 * <li>dscal_j
65 * <li>dswap_j
66 * <li>dnrm2_j
67 * <li>dcopy_j
68 * <li>drotg_j
69 * </ol>
70 * It also contains utility routines that the translator found useful
71 * while translating the FORTRAN code to Java code. "col" indicates that
72 * the routine operates on two columns of a matrix. "colv" indicates that
73 * the routine operates on a column of a matrix and a vector. The "p"
74 * at the end of dscalp, dnrm2p, and dcopyp indicates that these
75 * routines operate on a portion of a vector:
76 *
77 * <ol>
78 * <li>colisamax_j
79 * <li>colaxpy_j
80 * <li>colvaxpy_j
81 * <li>colvraxpy_j
82 * <li>coldot_j
83 * <li>colvdot_j
84 * <li>colscal_j
85 * <li>dscalp_j
86 * <li>colswap_j
87 * <li>colnrm2_j
88 * <li>dnrm2p_j
89 * <li>dcopyp_j
90 * <li>colrot_j
91 * <li>sign_j
92 * </ol>
93 *
94 * <p>
95 * <b>IMPORTANT:</b> The "_j" suffixes indicate that these routines use
96 * Java style indexing. For example, you will see
97 *
98 * for (i = 0; i < n; i++)
99 *
100 * rather than (FORTRAN style)
101 *
102 * for (i = 1; i <= n; i++)
103 *
104 * To use the "_j" routines you will have to
105 * fill elements 0 through n - 1 of vectors rather than elements 1
106 * through n.
107 * [Also, before using the isamax and colisamax methods
108 * make sure that they are doing what you expect them to do.]
109 * Versions of these programs that use FORTRAN style indexing are
110 * also available. They end with the suffix "_f77".
111 *
112 * <p>
113 * This class was translated by a statistician from FORTRAN versions of
114 * the LINPACK blas. It is NOT an official translation. When
115 * public domain Java numerical analysis routines become available
116 * from the people who produce LAPACK, then <b>THE CODE PRODUCED
117 * BY THE NUMERICAL ANALYSTS SHOULD BE USED</b>.
118 *
119 * <p>
120 * Meanwhile, if you have suggestions for improving this
121 * code, please contact Steve Verrill at sverrill@fs.fed.us.
122 *
123 * @author Steve Verrill
124 * @version .5 --- June 3, 1997
125 *
126 */
127 public class Blas extends Object {
128
129 /**
130 *
131 * <p>
132 * This method multiplies a constant times a portion of a column
133 * of a matrix and adds the product to the corresponding portion
134 * of another column of the matrix --- a portion of col2 is
135 * replaced by the corresponding portion of a*col1 + col2.
136 * It uses unrolled loops.
137 * It is a modification of the LINPACK subroutine
138 * DAXPY. In the LINPACK listing DAXPY is attributed to Jack Dongarra
139 * with a date of 3/11/78.
140 *
141 * Translated and modified by Steve Verrill, February 26, 1997.
142 *
143 * @param nrow The number of rows involved
144 * @param a The constant
145 * @param x[ ][ ] The matrix
146 * @param begin The starting row
147 * @param j1 The id of col1
148 * @param j2 The id of col2
149 *
150 */
151 public static void colaxpy_j( int nrow, double a, double x[][], int begin,
152 int j1, int j2 ) {
153
154 int i, m, mpbegin, end;
155
156 if ( nrow <= 0 ) return;
157 if ( a == 0.0 ) return;
158
159 m = nrow % 4;
160 mpbegin = m + begin;
161 end = begin + nrow - 1;
162
163 for ( i = begin; i < mpbegin; i++ ) {
164
165 x[i][j2] += a * x[i][j1];
166
167 }
168
169 for ( i = mpbegin; i <= end; i += 4 ) {
170
171 x[i][j2] += a * x[i][j1];
172 x[i + 1][j2] += a * x[i + 1][j1];
173 x[i + 2][j2] += a * x[i + 2][j1];
174 x[i + 3][j2] += a * x[i + 3][j1];
175
176 }
177
178 return;
179
180 }
181
182 /**
183 *
184 * <p>
185 * This method calculates the dot product of portions of two
186 * columns of a matrix. It uses unrolled loops.
187 * It is a modification of the LINPACK function
188 * DDOT. In the LINPACK listing DDOT is attributed to Jack Dongarra
189 * with a date of 3/11/78.
190 *
191 * Translated and modified by Steve Verrill, February 27, 1997.
192 *
193 * @param nrow The number of rows involved
194 * @param x[ ][ ] The matrix
195 * @param begin The starting row
196 * @param j1 The id of the first column
197 * @param j2 The id of the second column
198 *
199 */
200 public static double coldot_j( int nrow, double x[][], int begin,
201 int j1, int j2 ) {
202
203 double coldot;
204 int i, m, mpbegin, end;
205
206 coldot = 0.0;
207
208 if ( nrow <= 0 ) return coldot;
209
210 m = nrow % 5;
211 mpbegin = m + begin;
212 end = begin + nrow - 1;
213
214 for ( i = begin; i < mpbegin; i++ ) {
215
216 coldot += x[i][j1] * x[i][j2];
217
218 }
219
220 for ( i = mpbegin; i <= end; i += 5 ) {
221
222 coldot += x[i][j1] * x[i][j2] +
223 x[i + 1][j1] * x[i + 1][j2] +
224 x[i + 2][j1] * x[i + 2][j2] +
225 x[i + 3][j1] * x[i + 3][j2] +
226 x[i + 4][j1] * x[i + 4][j2];
227
228 }
229
230 return coldot;
231
232 }
233
234 /**
235 *
236 * <p>
237 * This method finds the index of the element of a portion of a
238 * column of a matrix that has the maximum absolute value.
239 * It is a modification of the LINPACK function ISAMAX.
240 * In the LINPACK listing ISAMAX is attributed to Jack Dongarra
241 * with a date of March 11, 1978.
242 *
243 * Translated by Steve Verrill, March 10, 1998.
244 *
245 * @param n The number of elements to be checked
246 * @param x[ ][ ] The matrix
247 * @param incx The subscript increment for x[ ][ ]
248 * @param begin The starting row
249 * @param j The id of the column
250 *
251 */
252 public static int colisamax_j( int n, double x[][], int incx,
253 int begin, int j ) {
254 //NOTE THAT THERE IS NO DIFFERENCE BETWEEN
255 //colisamax_j AND colisamax_f77. THIS IS NOT A MISTAKE.
256 double xmax;
257 int isamax, i, ix;
258
259 if ( n < 1 ) {
260
261 isamax = 0;
262
263 } else if ( n == 1 ) {
264
265 isamax = 1;
266
267 } else if ( incx == 1 ) {
268
269 isamax = 1;
270 ix = begin;
271 xmax = Math.abs( x[ix][j] );
272 ix++;
273
274 for ( i = 2; i <= n; i++ ) {
275
276 if ( Math.abs( x[ix][j] ) > xmax ) {
277
278 isamax = i;
279 xmax = Math.abs( x[ix][j] );
280
281 }
282
283 ix++;
284
285 }
286
287 } else {
288
289 isamax = 1;
290 ix = begin;
291 xmax = Math.abs( x[ix][j] );
292 ix += incx;
293
294 for ( i = 2; i <= n; i++ ) {
295
296 if ( Math.abs( x[ix][j] ) > xmax ) {
297
298 isamax = i;
299 xmax = Math.abs( x[ix][j] );
300
301 }
302
303 ix += incx;
304
305 }
306
307 }
308
309 return isamax;
310
311 }
312
313 /**
314 *
315 * <p>
316 * This method calculates the Euclidean norm of a portion of a
317 * column of a matrix.
318 * It is a modification of the LINPACK function
319 * dnrm2.
320 * In the LINPACK listing dnrm2 is attributed to C.L. Lawson
321 * with a date of January 8, 1978. The routine below is based
322 * on a more recent dnrm2 version that is attributed in LAPACK
323 * documentation to Sven Hammarling.
324 *
325 * Translated and modified by Steve Verrill, February 26, 1997.
326 *
327 * @param nrow The number of rows involved
328 * @param x[ ][ ] The matrix
329 * @param begin The starting row
330 * @param j The id of the column
331 *
332 */
333 public static double colnrm2_j( int nrow, double x[][], int begin, int j ) {
334
335 double absxij, norm, scale, ssq, fac;
336 int i, end;
337
338 if ( nrow < 1 ) {
339
340 norm = 0.0;
341
342 } else if ( nrow == 1 ) {
343
344 norm = Math.abs( x[begin][j] );
345
346 } else {
347
348 scale = 0.0;
349 ssq = 1.0;
350
351 end = begin + nrow - 1;
352
353 for ( i = begin; i <= end; i++ ) {
354
355 if ( x[i][j] != 0.0 ) {
356
357 absxij = Math.abs( x[i][j] );
358
359 if ( scale < absxij ) {
360
361 fac = scale / absxij;
362 ssq = 1.0 + ssq * fac * fac;
363 scale = absxij;
364
365 } else {
366
367 fac = absxij / scale;
368 ssq += fac * fac;
369
370 }
371
372 }
373
374 }
375
376 norm = scale * Math.sqrt( ssq );
377
378 }
379
380 return norm;
381
382 }
383
384 /**
385 *
386 * <p>
387 * This method "applies a plane rotation."
388 * It is a modification of the LINPACK function
389 * DROT. In the LINPACK listing DROT is attributed to Jack Dongarra
390 * with a date of 3/11/78.
391 *
392 * Translated and modified by Steve Verrill, March 4, 1997.
393 *
394 * @param n The order of x[ ][ ]
395 * @param x[ ][ ] The matrix
396 * @param j1 The id of the first column
397 * @param j2 The id of the second column
398 * @param c "cos"
399 * @param s "sin"
400 *
401 */
402
403 public static void colrot_j( int n, double x[][],
404 int j1, int j2, double c, double s ) {
405
406 double temp;
407 int i;
408
409 if ( n <= 0 ) return;
410
411 for ( i = 0; i < n; i++ ) {
412
413 temp = c * x[i][j1] + s * x[i][j2];
414 x[i][j2] = c * x[i][j2] - s * x[i][j1];
415 x[i][j1] = temp;
416
417 }
418
419 return;
420
421 }
422
423 /**
424 *
425 * <p>
426 * This method scales a portion of a column of a matrix by a constant.
427 * It uses unrolled loops.
428 * It is a modification of the LINPACK subroutine
429 * DSCAL. In the LINPACK listing DSCAL is attributed to Jack Dongarra
430 * with a date of 3/11/78.
431 *
432 * Translated and modified by Steve Verrill, February 27, 1997.
433 *
434 * @param nrow The number of rows involved
435 * @param a The constant
436 * @param x[ ][ ] The matrix
437 * @param begin The starting row
438 * @param j The id of the column
439 *
440 */
441
442 public static void colscal_j( int nrow, double a, double x[][], int begin, int j ) {
443
444 int i, m, mpbegin, end;
445
446 if ( nrow <= 0 ) return;
447
448 m = nrow % 5;
449 mpbegin = m + begin;
450 end = begin + nrow - 1;
451
452 for ( i = begin; i < mpbegin; i++ ) {
453
454 x[i][j] *= a;
455
456 }
457
458 for ( i = mpbegin; i <= end; i += 5 ) {
459
460 x[i][j] *= a;
461 x[i + 1][j] *= a;
462 x[i + 2][j] *= a;
463 x[i + 3][j] *= a;
464 x[i + 4][j] *= a;
465
466 }
467
468 return;
469
470 }
471
472 /**
473 *
474 * <p>
475 * This method interchanges two columns of a matrix.
476 * It uses unrolled loops.
477 * It is a modification of the LINPACK function
478 * DSWAP. In the LINPACK listing DSWAP is attributed to Jack Dongarra
479 * with a date of 3/11/78.
480 *
481 * Translated and modified by Steve Verrill, February 26, 1997.
482 *
483 * @param n The number of rows of the matrix
484 * @param x[ ][ ] The matrix
485 * @param j1 The id of the first column
486 * @param j2 The id of the second column
487 *
488 */
489 public static void colswap_j( int n, double x[][], int j1, int j2 ) {
490
491 double temp;
492 int i, m;
493
494 if ( n <= 0 ) return;
495
496 m = n % 3;
497
498 for ( i = 0; i < m; i++ ) {
499
500 temp = x[i][j1];
501 x[i][j1] = x[i][j2];
502 x[i][j2] = temp;
503
504 }
505
506 for ( i = m; i < n; i += 3 ) {
507
508 temp = x[i][j1];
509 x[i][j1] = x[i][j2];
510 x[i][j2] = temp;
511
512 temp = x[i + 1][j1];
513 x[i + 1][j1] = x[i + 1][j2];
514 x[i + 1][j2] = temp;
515
516 temp = x[i + 2][j1];
517 x[i + 2][j1] = x[i + 2][j2];
518 x[i + 2][j2] = temp;
519
520 }
521
522 return;
523
524 }
525
526 /**
527 *
528 * <p>
529 * This method multiplies a constant times a portion of a column
530 * of a matrix x[ ][ ] and adds the product to the corresponding portion
531 * of a vector y[ ] --- a portion of y[ ] is replaced by the corresponding
532 * portion of ax[ ][j] + y[ ].
533 * It uses unrolled loops.
534 * It is a modification of the LINPACK subroutine
535 * DAXPY. In the LINPACK listing DAXPY is attributed to Jack Dongarra
536 * with a date of 3/11/78.
537 *
538 * Translated and modified by Steve Verrill, March 1, 1997.
539 *
540 * @param nrow The number of rows involved
541 * @param a The constant
542 * @param x[ ][ ] The matrix
543 * @param y[ ] The vector
544 * @param begin The starting row
545 * @param j The id of the column of the x matrix
546 *
547 */
548
549 public static void colvaxpy_j( int nrow, double a, double x[][], double y[],
550 int begin, int j ) {
551
552 int i, m, mpbegin, end;
553
554 if ( nrow <= 0 ) return;
555 if ( a == 0.0 ) return;
556
557 m = nrow % 4;
558 mpbegin = m + begin;
559 end = begin + nrow - 1;
560
561 for ( i = begin; i < mpbegin; i++ ) {
562
563 y[i] += a * x[i][j];
564
565 }
566
567 for ( i = mpbegin; i <= end; i += 4 ) {
568
569 y[i] += a * x[i][j];
570 y[i + 1] += a * x[i + 1][j];
571 y[i + 2] += a * x[i + 2][j];
572 y[i + 3] += a * x[i + 3][j];
573
574 }
575
576 return;
577
578 }
579
580 /**
581 *
582 * <p>
583 * This method calculates the dot product of a portion of a
584 * column of a matrix and the corresponding portion of a
585 * vector. It uses unrolled loops.
586 * It is a modification of the LINPACK function
587 * DDOT. In the LINPACK listing DDOT is attributed to Jack Dongarra
588 * with a date of 3/11/78.
589 *
590 * Translated and modified by Steve Verrill, March 1, 1997.
591 *
592 * @param nrow The number of rows involved
593 * @param x[ ][ ] The matrix
594 * @param y[ ] The vector
595 * @param begin The starting row
596 * @param j The id of the column of the matrix
597 *
598 */
599
600 public static double colvdot_j( int nrow, double x[][], double y[],
601 int begin, int j ) {
602
603 double colvdot;
604 int i, m, mpbegin, end;
605
606 colvdot = 0.0;
607
608 if ( nrow <= 0 ) return colvdot;
609
610 m = nrow % 5;
611 mpbegin = m + begin;
612 end = begin + nrow - 1;
613
614 for ( i = begin; i < mpbegin; i++ ) {
615
616 colvdot += x[i][j] * y[i];
617
618 }
619
620 for ( i = mpbegin; i <= end; i += 5 ) {
621
622 colvdot += x[i][j] * y[i] +
623 x[i + 1][j] * y[i + 1] +
624 x[i + 2][j] * y[i + 2] +
625 x[i + 3][j] * y[i + 3] +
626 x[i + 4][j] * y[i + 4];
627
628 }
629
630 return colvdot;
631
632 }
633
634 /**
635 *
636 * <p>
637 * This method multiplies a constant times a portion of a vector y[ ]
638 * and adds the product to the corresponding portion
639 * of a column of a matrix x[ ][ ] --- a portion of column j of x[ ][ ]
640 * is replaced by the corresponding
641 * portion of ay[ ] + x[ ][j].
642 * It uses unrolled loops.
643 * It is a modification of the LINPACK subroutine
644 * DAXPY. In the LINPACK listing DAXPY is attributed to Jack Dongarra
645 * with a date of 3/11/78.
646 *
647 * Translated and modified by Steve Verrill, March 3, 1997.
648 *
649 * @param nrow The number of rows involved
650 * @param a The constant
651 * @param y[ ] The vector
652 * @param x[ ][ ] The matrix
653 * @param begin The starting row
654 * @param j The id of the column of the x matrix
655 *
656 */
657
658 public static void colvraxpy_j( int nrow, double a, double y[], double x[][],
659 int begin, int j ) {
660
661 int i, m, mpbegin, end;
662
663 if ( nrow <= 0 ) return;
664 if ( a == 0.0 ) return;
665
666 m = nrow % 4;
667 mpbegin = m + begin;
668 end = begin + nrow - 1;
669
670 for ( i = begin; i < mpbegin; i++ ) {
671
672 x[i][j] += a * y[i];
673
674 }
675
676 for ( i = mpbegin; i <= end; i += 4 ) {
677
678 x[i][j] += a * y[i];
679 x[i + 1][j] += a * y[i + 1];
680 x[i + 2][j] += a * y[i + 2];
681 x[i + 3][j] += a * y[i + 3];
682
683 }
684
685 return;
686
687 }
688
689 /**
690 *
691 * <p>
692 * This method multiplies a constant times a vector and adds the product
693 * to another vector --- dy[ ] = da*dx[ ] + dy[ ].
694 * It uses unrolled loops for increments equal to
695 * one. It is a translation from FORTRAN to Java of the LINPACK subroutine
696 * DAXPY. In the LINPACK listing DAXPY is attributed to Jack Dongarra
697 * with a date of 3/11/78.
698 *
699 * Translated by Steve Verrill, June 3, 1997.
700 *
701 * @param n The order of the vectors dy[ ] and dx[ ]
702 * @param da The constant
703 * @param dx[ ] This vector will be multiplied by the constant da
704 * @param incx The subscript increment for dx[ ]
705 * @param dy[ ] This vector will be added to da*dx[ ]
706 * @param incy The subscript increment for dy[ ]
707 *
708 */
709
710 public static void daxpy_j( int n, double da, double dx[], int incx, double dy[], int incy ) {
711
712 int i, ix, iy, m;
713
714 if ( n <= 0 ) return;
715 if ( da == 0.0 ) return;
716
717 if ( ( incx == 1 ) && ( incy == 1 ) ) {
718
719 //both increments equal to 1
720
721 m = n % 4;
722
723 for ( i = 0; i < m; i++ ) {
724
725 dy[i] += da * dx[i];
726
727 }
728
729 for ( i = m; i < n; i += 4 ) {
730
731 dy[i] += da * dx[i];
732 dy[i + 1] += da * dx[i + 1];
733 dy[i + 2] += da * dx[i + 2];
734 dy[i + 3] += da * dx[i + 3];
735
736 }
737
738 return;
739
740 }
741
742 //at least one increment not equal to 1
743
744 ix = 0;
745 iy = 0;
746
747 if ( incx < 0 ) ix = ( -n + 1 ) * incx;
748 if ( incy < 0 ) iy = ( -n + 1 ) * incy;
749
750 for ( i = 0; i < n; i++ ) {
751
752 dy[iy] += da * dx[ix];
753
754 ix += incx;
755 iy += incy;
756
757 }
758
759 return;
760
761 }
762
763 /**
764 *
765 * <p>
766 * This method copies the vector dx[ ] to the vector dy[ ].
767 * It uses unrolled loops for increments equal to
768 * one. It is a translation from FORTRAN to Java of the LINPACK subroutine
769 * DCOPY. In the LINPACK listing DCOPY is attributed to Jack Dongarra
770 * with a date of 3/11/78.
771 *
772 * Translated by Steve Verrill, March 1, 1997.
773 *
774 * @param n The order of dx[ ] and dy[ ]
775 * @param dx[ ] vector
776 * @param incx The subscript increment for dx[ ]
777 * @param dy[ ] vector
778 * @param incy The subscript increment for dy[ ]
779 *
780 */
781 public static void dcopy_j( int n, double dx[], int incx, double dy[], int incy ) {
782
783 // double dtemp;
784 int i, ix, iy, m;
785
786 if ( n <= 0 ) return;
787
788 if ( ( incx == 1 ) && ( incy == 1 ) ) {
789
790 //both increments equal to 1
791
792 m = n % 7;
793
794 for ( i = 0; i < m; i++ ) {
795
796 dy[i] = dx[i];
797
798 }
799
800 for ( i = m; i < n; i += 7 ) {
801
802 dy[i] = dx[i];
803 dy[i + 1] = dx[i + 1];
804 dy[i + 2] = dx[i + 2];
805 dy[i + 3] = dx[i + 3];
806 dy[i + 4] = dx[i + 4];
807 dy[i + 5] = dx[i + 5];
808 dy[i + 6] = dx[i + 6];
809
810 }
811
812 return;
813
814 }
815
816 //at least one increment not equal to 1
817
818 ix = 0;
819 iy = 0;
820
821 if ( incx < 0 ) ix = ( -n + 1 ) * incx;
822 if ( incy < 0 ) iy = ( -n + 1 ) * incy;
823
824 for ( i = 0; i < n; i++ ) {
825
826 dy[iy] = dx[ix];
827
828 ix += incx;
829 iy += incy;
830
831 }
832
833 return;
834
835 }
836
837 /**
838 *
839 * <p>
840 * This method copies a portion of vector x[ ] to the corresponding
841 * portion of vector y[ ].
842 * It uses unrolled loops.
843 * It is a modification of the LINPACK subroutine
844 * dcopy. In the LINPACK listing dcopy is attributed to Jack Dongarra
845 * with a date of 3/11/78.
846 *
847 * Translated by Steve Verrill, March 1, 1997.
848 *
849 * @param nrow The number of rows involved
850 * @param x[ ] vector
851 * @param y[ ] vector
852 * @param begin The starting row
853 *
854 */
855
856 public static void dcopyp_j( int nrow, double x[], double y[], int begin ) {
857
858 // double temp;
859 int i, m, mpbegin, end;
860
861 m = nrow % 7;
862 mpbegin = m + begin;
863 end = begin + nrow - 1;
864
865 for ( i = begin; i < mpbegin; i++ ) {
866
867 y[i] = x[i];
868
869 }
870
871 for ( i = mpbegin; i <= end; i += 7 ) {
872
873 y[i] = x[i];
874 y[i + 1] = x[i + 1];
875 y[i + 2] = x[i + 2];
876 y[i + 3] = x[i + 3];
877 y[i + 4] = x[i + 4];
878 y[i + 5] = x[i + 5];
879 y[i + 6] = x[i + 6];
880
881 }
882
883 return;
884
885 }
886
887 /**
888 *
889 * <p>
890 * This method calculates the dot product of two vectors.
891 * It uses unrolled loops for increments equal to
892 * one. It is a translation from FORTRAN to Java of the LINPACK function
893 * DDOT. In the LINPACK listing DDOT is attributed to Jack Dongarra
894 * with a date of 3/11/78.
895 *
896 * Translated by Steve Verrill, June 3, 1997.
897 *
898 * @param n The order of the vectors dx[ ] and dy[ ]
899 * @param dx[ ] vector
900 * @param incx The subscript increment for dx[ ]
901 * @param dy[ ] vector
902 * @param incy The subscript increment for dy[ ]
903 *
904 */
905 public static double ddot_j( int n, double dx[], int incx, double dy[], int incy ) {
906
907 double ddot;
908 int i, ix, iy, m;
909
910 ddot = 0.0;
911
912 if ( n <= 0 ) return ddot;
913
914 if ( ( incx == 1 ) && ( incy == 1 ) ) {
915
916 //both increments equal to 1
917
918 m = n % 5;
919
920 for ( i = 0; i < m; i++ ) {
921
922 ddot += dx[i] * dy[i];
923
924 }
925
926 for ( i = m; i < n; i += 5 ) {
927
928 ddot += dx[i] * dy[i] + dx[i + 1] * dy[i + 1] + dx[i + 2] * dy[i + 2] +
929 dx[i + 3] * dy[i + 3] + dx[i + 4] * dy[i + 4];
930
931 }
932
933 return ddot;
934
935 }
936
937 //at least one increment not equal to 1
938
939 ix = 0;
940 iy = 0;
941
942 if ( incx < 0 ) ix = ( -n + 1 ) * incx;
943 if ( incy < 0 ) iy = ( -n + 1 ) * incy;
944
945 for ( i = 0; i < n; i++ ) {
946
947 ddot += dx[ix] * dy[iy];
948
949 ix += incx;
950 iy += incy;
951
952 }
953
954 return ddot;
955
956 }
957
958 /**
959 *
960 * <p>
961 * This method calculates the Euclidean norm of the vector
962 * stored in dx[ ] with storage increment incx.
963 * It is a translation from FORTRAN to Java of the LINPACK function
964 * DNRM2.
965 * In the LINPACK listing DNRM2 is attributed to C.L. Lawson
966 * with a date of January 8, 1978. The routine below is based
967 * on a more recent DNRM2 version that is attributed in LAPACK
968 * documentation to Sven Hammarling.
969 *
970 * Translated by Steve Verrill, June 3, 1997.
971 *
972 * @param n The order of the vector x[ ]
973 * @param x[ ] vector
974 * @param incx The subscript increment for x[ ]
975 *
976 */
977
978 public static double dnrm2_j( int n, double x[], int incx ) {
979
980 double absxi, norm, scale, ssq, fac;
981 int ix, limit;
982
983 if ( n < 1 || incx < 1 ) {
984
985 norm = 0.0;
986
987 } else if ( n == 1 ) {
988
989 norm = Math.abs( x[0] );
990
991 } else {
992
993 scale = 0.0;
994 ssq = 1.0;
995
996 limit = ( n - 1 ) * incx;
997
998 for ( ix = 0; ix <= limit; ix += incx ) {
999
1000 if ( x[ix] != 0.0 ) {
1001
1002 absxi = Math.abs( x[ix] );
1003
1004 if ( scale < absxi ) {
1005
1006 fac = scale / absxi;
1007 ssq = 1.0 + ssq * fac * fac;
1008 scale = absxi;
1009
1010 } else {
1011
1012 fac = absxi / scale;
1013 ssq += fac * fac;
1014
1015 }
1016
1017 }
1018
1019 }
1020
1021 norm = scale * Math.sqrt( ssq );
1022
1023 }
1024
1025 return norm;
1026
1027 }
1028
1029 /**
1030 *
1031 * <p>
1032 * This method calculates the Euclidean norm of a portion
1033 * of a vector x[ ].
1034 * It is a modification of the LINPACK function
1035 * dnrm2.
1036 * In the LINPACK listing dnrm2 is attributed to C.L. Lawson
1037 * with a date of January 8, 1978. The routine below is based
1038 * on a more recent dnrm2 version that is attributed in LAPACK
1039 * documentation to Sven Hammarling.
1040 *
1041 * Translated by Steve Verrill, March 3, 1997.
1042 *
1043 * @param nrow The number of rows involved
1044 * @param x[ ] vector
1045 * @param begin The starting row
1046 *
1047 */
1048 public static double dnrm2p_j( int nrow, double x[], int begin ) {
1049
1050 double absxi, norm, scale, ssq, fac;
1051 int i, end;
1052
1053 if ( nrow < 1 ) {
1054
1055 norm = 0.0;
1056
1057 } else if ( nrow == 1 ) {
1058
1059 norm = Math.abs( x[begin] );
1060
1061 } else {
1062
1063 scale = 0.0;
1064 ssq = 1.0;
1065
1066 end = begin + nrow - 1;
1067
1068 for ( i = begin; i <= end; i++ ) {
1069
1070 if ( x[i] != 0.0 ) {
1071
1072 absxi = Math.abs( x[i] );
1073
1074 if ( scale < absxi ) {
1075
1076 fac = scale / absxi;
1077 ssq = 1.0 + ssq * fac * fac;
1078 scale = absxi;
1079
1080 } else {
1081
1082 fac = absxi / scale;
1083 ssq += fac * fac;
1084
1085 }
1086
1087 }
1088
1089 }
1090
1091 norm = scale * Math.sqrt( ssq );
1092
1093 }
1094
1095 return norm;
1096
1097 }
1098
1099 /**
1100 *
1101 * <p>
1102 * This method constructs a Givens plane rotation.
1103 * It is a translation from FORTRAN to Java of the LINPACK subroutine
1104 * DROTG. In the LINPACK listing DROTG is attributed to Jack Dongarra
1105 * with a date of 3/11/78.
1106 *
1107 * Translated by Steve Verrill, March 3, 1997.
1108 *
1109 * @param rotvec[] Contains the a,b,c,s values. In Java they
1110 * cannot be passed as primitive types (e.g., double
1111 * or int or ...) if we want their return values to
1112 * be altered.
1113 *
1114 */
1115 public static void drotg_j( double rotvec[] ) {
1116
1117 //construct a Givens plane rotation
1118
1119 double a, b, c, s, roe, scale, r, z, ra, rb;
1120
1121 a = rotvec[0];
1122 b = rotvec[1];
1123
1124 roe = b;
1125
1126 if ( Math.abs( a ) > Math.abs( b ) ) roe = a;
1127
1128 scale = Math.abs( a ) + Math.abs( b );
1129
1130 if ( scale != 0.0 ) {
1131
1132 ra = a / scale;
1133 rb = b / scale;
1134 r = scale * Math.sqrt( ra * ra + rb * rb );
1135 r = sign_j( 1.0, roe ) * r;
1136 c = a / r;
1137 s = b / r;
1138 z = 1.0;
1139 if ( Math.abs( a ) > Math.abs( b ) ) z = s;
1140 if ( ( Math.abs( b ) >= Math.abs( a ) ) && ( c != 0.0 ) ) z = 1.0 / c;
1141
1142 } else {
1143
1144 c = 1.0;
1145 s = 0.0;
1146 r = 0.0;
1147 z = 0.0;
1148
1149 }
1150
1151 a = r;
1152 b = z;
1153
1154 rotvec[0] = a;
1155 rotvec[1] = b;
1156 rotvec[2] = c;
1157 rotvec[3] = s;
1158
1159 return;
1160
1161 }
1162
1163 /**
1164 *
1165 * <p>
1166 * This method scales a vector by a constant.
1167 * It uses unrolled loops for an increment equal to
1168 * one. It is a translation from FORTRAN to Java of the LINPACK subroutine
1169 * DSCAL. In the LINPACK listing DSCAL is attributed to Jack Dongarra
1170 * with a date of 3/11/78.
1171 *
1172 * Translated by Steve Verrill, June 3, 1997.
1173 *
1174 * @param n The order of the vector dx[ ]
1175 * @param da The constant
1176 * @param dx[ ] This vector will be multiplied by the constant da
1177 * @param incx The subscript increment for dx[ ]
1178 *
1179 */
1180 public static void dscal_j( int n, double da, double dx[], int incx ) {
1181
1182 int i, m, nincx;
1183
1184 if ( n <= 0 || incx <= 0 ) return;
1185
1186 if ( incx == 1 ) {
1187
1188 //increment equal to 1
1189
1190 m = n % 5;
1191
1192 for ( i = 0; i < m; i++ ) {
1193
1194 dx[i] *= da;
1195
1196 }
1197
1198 for ( i = m; i < n; i += 5 ) {
1199
1200 dx[i] *= da;
1201 dx[i + 1] *= da;
1202 dx[i + 2] *= da;
1203 dx[i + 3] *= da;
1204 dx[i + 4] *= da;
1205
1206 }
1207
1208 return;
1209
1210 }
1211
1212 //increment not equal to 1
1213
1214 nincx = n * incx;
1215
1216 for ( i = 0; i < nincx; i += incx ) {
1217
1218 dx[i] *= da;
1219
1220 }
1221
1222 return;
1223
1224 }
1225
1226 /**
1227 *
1228 * <p>
1229 * This method scales a portion of a vector by a constant.
1230 * It uses unrolled loops.
1231 * It is a modification of the LINPACK subroutine
1232 * DSCAL. In the LINPACK listing DSCAL is attributed to Jack Dongarra
1233 * with a date of 3/11/78.
1234 *
1235 * Translated and modified by Steve Verrill, March 3, 1997.
1236 *
1237 * @param nrow The number of rows involved
1238 * @param a The constant
1239 * @param x[ ] The vector
1240 * @param begin The starting row
1241 *
1242 */
1243
1244 public static void dscalp_j( int nrow, double a, double x[], int begin ) {
1245
1246 int i, m, mpbegin, end;
1247
1248 if ( nrow <= 0 ) return;
1249
1250 m = nrow % 5;
1251 mpbegin = m + begin;
1252 end = begin + nrow - 1;
1253
1254 for ( i = begin; i < mpbegin; i++ ) {
1255
1256 x[i] *= a;
1257
1258 }
1259
1260 for ( i = mpbegin; i <= end; i += 5 ) {
1261
1262 x[i] *= a;
1263 x[i + 1] *= a;
1264 x[i + 2] *= a;
1265 x[i + 3] *= a;
1266 x[i + 4] *= a;
1267
1268 }
1269
1270 return;
1271
1272 }
1273
1274 /**
1275 *
1276 * <p>
1277 * This method interchanges two vectors.
1278 * It uses unrolled loops for increments equal to
1279 * one. It is a translation from FORTRAN to Java of the LINPACK function
1280 * DSWAP. In the LINPACK listing DSWAP is attributed to Jack Dongarra
1281 * with a date of 3/11/78.
1282 *
1283 * Translated by Steve Verrill, June 3, 1997.
1284 *
1285 * @param n The order of the vectors dx[ ] and dy[ ]
1286 * @param dx[ ] vector
1287 * @param incx The subscript increment for dx[ ]
1288 * @param dy[ ] vector
1289 * @param incy The subscript increment for dy[ ]
1290 *
1291 */
1292
1293 public static void dswap_j( int n, double dx[], int incx, double dy[], int incy ) {
1294
1295 double dtemp;
1296 int i, ix, iy, m;
1297
1298 if ( n <= 0 ) return;
1299
1300 if ( ( incx == 1 ) && ( incy == 1 ) ) {
1301
1302 //both increments equal to 1
1303
1304 m = n % 3;
1305
1306 for ( i = 0; i < m; i++ ) {
1307
1308 dtemp = dx[i];
1309 dx[i] = dy[i];
1310 dy[i] = dtemp;
1311
1312 }
1313
1314 for ( i = m; i < n; i += 3 ) {
1315
1316 dtemp = dx[i];
1317 dx[i] = dy[i];
1318 dy[i] = dtemp;
1319
1320 dtemp = dx[i + 1];
1321 dx[i + 1] = dy[i + 1];
1322 dy[i + 1] = dtemp;
1323
1324 dtemp = dx[i + 2];
1325 dx[i + 2] = dy[i + 2];
1326 dy[i + 2] = dtemp;
1327
1328 }
1329
1330 return;
1331
1332 }
1333
1334 //at least one increment not equal to 1
1335
1336 ix = 0;
1337 iy = 0;
1338
1339 if ( incx < 0 ) ix = ( -n + 1 ) * incx;
1340 if ( incy < 0 ) iy = ( -n + 1 ) * incy;
1341
1342 for ( i = 0; i < n; i++ ) {
1343
1344 dtemp = dx[ix];
1345 dx[ix] = dy[iy];
1346 dy[iy] = dtemp;
1347
1348 ix += incx;
1349 iy += incy;
1350
1351 }
1352
1353 return;
1354
1355 }
1356
1357 /**
1358 *
1359 * <p>
1360 * This method finds the index of the element of a vector
1361 * that has the maximum absolute value. It is a translation
1362 * from FORTRAN to Java of the LINPACK function ISAMAX.
1363 * In the LINPACK listing ISAMAX is attributed to Jack Dongarra
1364 * with a date of March 11, 1978. <b>Before you use this
1365 * version of isamax, make certain that it is doing what you
1366 * expect it to do.</b>
1367 *
1368 * Translated by Steve Verrill, March 10, 1998.
1369 *
1370 * @param n The number of elements to be checked
1371 * @param x[ ] vector
1372 * @param incx The subscript increment for x[ ]
1373 *
1374 */
1375 public static int isamax_j( int n, double x[], int incx ) {
1376
1377 double xmax;
1378 int isamax, i, ix;
1379
1380 if ( n < 1 ) {
1381
1382 isamax = 0;
1383
1384 } else if ( n == 1 ) {
1385
1386 isamax = 1;
1387
1388 } else if ( incx == 1 ) {
1389
1390 isamax = 1;
1391 xmax = Math.abs( x[0] );
1392
1393 for ( i = 1; i < n; i++ ) {
1394
1395 if ( Math.abs( x[i] ) > xmax ) {
1396
1397 isamax = i + 1;
1398 xmax = Math.abs( x[i] );
1399
1400 }
1401
1402 }
1403
1404 } else {
1405
1406 isamax = 1;
1407 ix = 0;
1408 xmax = Math.abs( x[ix] );
1409 ix += incx;
1410
1411 //The i = 2; i <= n material is correct here
1412 //because we are not indexing x[] by i.
1413
1414 for ( i = 2; i <= n; i++ ) {
1415
1416 if ( Math.abs( x[ix] ) > xmax ) {
1417
1418 isamax = i;
1419 xmax = Math.abs( x[ix] );
1420
1421 }
1422
1423 ix += incx;
1424
1425 }
1426
1427 }
1428
1429 return isamax;
1430
1431 }
1432
1433 /**
1434 *
1435 * <p>
1436 * This method multiplies an n x p matrix by a p x r matrix.
1437 *
1438 * Created by Steve Verrill, March 1997.
1439 *
1440 * @param a[ ][ ] The left matrix
1441 * @param b[ ][ ] The right matrix
1442 * @param c[ ][ ] The product
1443 * @param n n
1444 * @param p p
1445 * @param r r
1446 *
1447 */
1448 public static void matmat_j( double a[][], double b[][], double c[][],
1449 int n, int p, int r ) {
1450
1451 double vdot;
1452 int i, j, k, m;
1453
1454 for ( i = 0; i < n; i++ ) {
1455
1456 for ( j = 0; j < r; j++ ) {
1457
1458 vdot = 0.0;
1459
1460 m = p % 5;
1461
1462 for ( k = 0; k < m; k++ ) {
1463
1464 vdot += a[i][k] * b[k][j];
1465
1466 }
1467
1468 for ( k = m; k < p; k += 5 ) {
1469
1470 vdot += a[i][k] * b[k][j] +
1471 a[i][k + 1] * b[k + 1][j] +
1472 a[i][k + 2] * b[k + 2][j] +
1473 a[i][k + 3] * b[k + 3][j] +
1474 a[i][k + 4] * b[k + 4][j];
1475
1476 }
1477
1478 c[i][j] = vdot;
1479
1480 }
1481
1482 }
1483
1484 }
1485
1486 /**
1487 *
1488 * This method obtains the transpose of an n x p matrix.
1489 *
1490 * Created by Steve Verrill, March 1997.
1491 *
1492 * @param a[ ][ ] matrix
1493 * @param at[ ][ ] transpose of the matrix
1494 * @param n n
1495 * @param p p
1496 *
1497 */
1498 public static void mattran_j( double a[][], double at[][],
1499 int n, int p ) {
1500
1501 int i, j;
1502
1503 for ( i = 0; i < n; i++ ) {
1504
1505 for ( j = 0; j < p; j++ ) {
1506
1507 at[j][i] = a[i][j];
1508
1509 }
1510
1511 }
1512
1513 }
1514
1515 /**
1516 *
1517 * This method multiplies an n x p matrix by a p x 1 vector.
1518 *
1519 * Created by Steve Verrill, March 1997.
1520 *
1521 * @param a[ ][ ] The matrix
1522 * @param b[ ] The vector
1523 * @param c[ ] The product
1524 * @param n n
1525 * @param p p
1526 *
1527 */
1528 public static void matvec_j( double a[][], double b[], double c[],
1529 int n, int p ) {
1530
1531 double vdot;
1532 int i, j, m;
1533
1534 for ( i = 0; i < n; i++ ) {
1535
1536 vdot = 0.0;
1537
1538 m = p % 5;
1539
1540 for ( j = 0; j < m; j++ ) {
1541
1542 vdot += a[i][j] * b[j];
1543
1544 }
1545
1546 for ( j = m; j < p; j += 5 ) {
1547
1548 vdot += a[i][j] * b[j] +
1549 a[i][j + 1] * b[j + 1] +
1550 a[i][j + 2] * b[j + 2] +
1551 a[i][j + 3] * b[j + 3] +
1552 a[i][j + 4] * b[j + 4];
1553
1554 }
1555
1556 c[i] = vdot;
1557
1558 }
1559
1560 }
1561
1562 /**
1563 *
1564 * This method implements the FORTRAN sign (not sin) function.
1565 * See the code for details.
1566 *
1567 * Created by Steve Verrill, March 1997.
1568 *
1569 * @param a a
1570 * @param b b
1571 *
1572 */
1573 public static double sign_j( double a, double b ) {
1574
1575 if ( b < 0.0 ) {
1576 return -Math.abs( a );
1577 }
1578
1579 return Math.abs( a );
1580
1581 }
1582
1583 }