1 ///////////////////////////////////////////////////////////////////////////////
2 // Copyright 2011 John Maddock. Distributed under the Boost
3 // Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 /* 1000d.f -- translated by f2c (version 20050501).
7 You must link the resulting object file with libf2c:
8 on Microsoft Windows system, link with libf2c.lib;
9 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
10 or, if you install libf2c.a in a standard place, with -lf2c -lm
11 -- in that order, at the end of the command line, as in
13 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
15 http://www.netlib.org/f2c/libf2c.zip
21 #if defined(TEST_GMPXX)
23 typedef mpf_class real_type
;
24 #elif defined(TEST_MPFRXX)
26 typedef mpfr_class real_type
;
27 #elif defined(TEST_CPP_DEC_FLOAT)
28 #include <boost/multiprecision/cpp_dec_float.hpp>
29 typedef boost::multiprecision::cpp_dec_float_50 real_type
;
30 #elif defined(TEST_MPFR_50)
31 #include <boost/multiprecision/mpfr.hpp>
32 typedef boost::multiprecision::mpfr_float_50 real_type
;
33 #elif defined(TEST_MPF_50)
34 #include <boost/multiprecision/gmp.hpp>
35 typedef boost::multiprecision::mpf_float_50 real_type
;
36 #elif defined(NATIVE_FLOAT128)
37 #include <boost/multiprecision/float128.hpp>
38 typedef __float128 real_type
;
40 std::ostream
& operator<<(std::ostream
& os
, const __float128
& f
)
42 return os
<< boost::multiprecision::float128(f
);
45 #include <boost/type_traits/has_left_shift.hpp>
50 struct has_left_shift
<std::basic_ostream
<char>, __float128
> : public mpl::true_
54 double lexical_cast
<double, __float128
>(const __float128
& f
)
61 #elif defined(TEST_FLOAT128)
62 #include <boost/multiprecision/float128.hpp>
63 typedef boost::multiprecision::float128 real_type
;
64 #elif defined(TEST_CPP_BIN_FLOAT_QUAD)
65 #include <boost/multiprecision/cpp_bin_float.hpp>
66 typedef boost::multiprecision::cpp_bin_float_quad real_type
;
67 #elif defined(TEST_CPP_BIN_FLOAT_OCT)
68 #include <boost/multiprecision/cpp_bin_float.hpp>
69 typedef boost::multiprecision::cpp_bin_float_oct real_type
;
71 typedef double real_type
;
74 #include <boost/lexical_cast.hpp>
77 #define CAST_TO_RT(x) x
82 integer
s_wsfe(cilist
*), e_wsfe(void), do_fio(integer
*, char*, ftnlen
),
83 s_wsle(cilist
*), do_lio(integer
*, integer
*, char*, ftnlen
),
85 /* Subroutine */ int s_stop(char*, ftnlen
);
100 /* Table of constant values */
102 static integer c__0
= 0;
103 static real_type c_b7
= CAST_TO_RT(1);
104 static integer c__1
= 1;
105 static integer c__9
= 9;
107 inline double second_(void)
109 return ((double)(clock())) / CLOCKS_PER_SEC
;
112 int dgefa_(real_type
*, integer
*, integer
*, integer
*, integer
*), dgesl_(real_type
*, integer
*, integer
*, integer
*, real_type
*, integer
*);
113 int dmxpy_(integer
*, real_type
*, integer
*, integer
*, real_type
*, real_type
*);
114 int matgen_(real_type
*, integer
*, integer
*, real_type
*, real_type
*);
115 real_type
epslon_(real_type
*);
116 real_type
ran_(integer
*);
117 int dscal_(integer
*, real_type
*, real_type
*, integer
*);
118 int daxpy_(integer
*, real_type
*, real_type
*, integer
*, real_type
*, integer
*);
119 integer
idamax_(integer
*, real_type
*, integer
*);
120 real_type
ddot_(integer
*, real_type
*, integer
*, real_type
*, integer
*);
121 int daxpy_(integer
*, real_type
*, real_type
*, integer
*, real_type
*, integer
*);
122 int dmxpy_(integer
*, real_type
*, integer
*, integer
*, real_type
*, real_type
*);
124 extern "C" int MAIN__()
127 std::cout
<< "Testing number<mpf_float<50> >" << std::endl
;
128 #elif defined(TEST_MPFR_50)
129 std::cout
<< "Testing number<mpf_float<50> >" << std::endl
;
130 #elif defined(TEST_GMPXX)
131 std::cout
<< "Testing mpf_class at 50 decimal degits" << std::endl
;
132 mpf_set_default_prec(((50 + 1) * 1000L) / 301L);
133 #elif defined(TEST_MPFRXX)
134 std::cout
<< "Testing mpfr_class at 50 decimal degits" << std::endl
;
135 mpfr_set_default_prec(((50 + 1) * 1000L) / 301L);
136 #elif defined(TEST_CPP_DEC_FLOAT)
137 std::cout
<< "Testing number<cpp_dec_float<50> >" << std::endl
;
138 #elif defined(NATIVE_FLOAT128)
139 std::cout
<< "Testing __float128" << std::endl
;
140 #elif defined(TEST_FLOAT128)
141 std::cout
<< "Testing number<float128_backend, et_off>" << std::endl
;
143 std::cout
<< "Testing double" << std::endl
;
147 static char fmt_1
[] = "(\002 Please send the results of this run to:\002"
148 "//\002 Jack J. Dongarra\002/\002 Computer Science Department\002/"
149 "\002 University of Tennessee\002/\002 Knoxville, Tennessee 37996"
150 "-1300\002//\002 Fax: 615-974-8296\002//\002 Internet: dongarra@c"
152 static char fmt_40
[] = "(\002 norm. resid resid mac"
153 "hep\002,\002 x(1) x(n)\002)";
154 static char fmt_50
[] = "(1p5e16.8)";
155 static char fmt_60
[] = "(//\002 times are reported for matrices of or"
157 static char fmt_70
[] = "(6x,\002factor\002,5x,\002solve\002,6x,\002tota"
158 "l\002,5x,\002mflops\002,7x,\002unit\002,6x,\002ratio\002)";
159 static char fmt_80
[] = "(\002 times for array with leading dimension o"
161 static char fmt_110
[] = "(6(1pe11.3))";
163 /* System generated locals */
165 real_type d__1
, d__2
, d__3
;
167 /* Builtin functions */
169 /* Local variables */
170 static real_type a
[1001000] /* was [1001][1000] */, b
[1000];
171 static integer i__
, n
;
172 static real_type x
[1000];
176 static real_type eps
;
178 static double time
[6], cray
, total
;
179 static integer ipvt
[1000];
180 static real_type resid
, norma
;
181 static real_type normx
;
182 static real_type residn
;
184 /* Fortran I/O blocks */
185 static cilist io___4
= {0, 6, 0, fmt_1
, 0};
186 static cilist io___20
= {0, 6, 0, fmt_40
, 0};
187 static cilist io___21
= {0, 6, 0, fmt_50
, 0};
188 static cilist io___22
= {0, 6, 0, fmt_60
, 0};
189 static cilist io___23
= {0, 6, 0, fmt_70
, 0};
190 static cilist io___24
= {0, 6, 0, fmt_80
, 0};
191 static cilist io___25
= {0, 6, 0, fmt_110
, 0};
192 static cilist io___26
= {0, 6, 0, 0, 0};
196 /* this program was updated on 10/12/92 to correct a */
197 /* problem with the random number generator. The previous */
198 /* random number generator had a short period and produced */
199 /* singular matrices occasionally. */
205 /* Computing 3rd power */
207 /* Computing 2nd power */
209 ops
= boost::lexical_cast
<double>(real_type(d__1
* (d__1
* d__1
) * 2. / 3. + d__2
* d__2
* 2.));
211 matgen_(a
, &lda
, &n
, b
, &norma
);
213 /* ****************************************************************** */
214 /* ****************************************************************** */
215 /* you should replace the call to dgefa and dgesl */
216 /* by calls to your linear equation solver. */
217 /* ****************************************************************** */
218 /* ****************************************************************** */
221 dgefa_(a
, &lda
, &n
, ipvt
, &info
);
222 time
[0] = second_() - t1
;
224 dgesl_(a
, &lda
, &n
, ipvt
, b
, &c__0
);
225 time
[1] = second_() - t1
;
226 total
= time
[0] + time
[1];
227 /* ****************************************************************** */
228 /* ****************************************************************** */
230 /* compute a residual to verify results. */
233 for (i__
= 1; i__
<= i__1
; ++i__
)
235 x
[i__
- 1] = b
[i__
- 1];
238 matgen_(a
, &lda
, &n
, b
, &norma
);
240 for (i__
= 1; i__
<= i__1
; ++i__
)
242 b
[i__
- 1] = -b
[i__
- 1];
245 dmxpy_(&n
, b
, &n
, &lda
, x
, a
);
246 resid
= CAST_TO_RT(0);
247 normx
= CAST_TO_RT(0);
249 for (i__
= 1; i__
<= i__1
; ++i__
)
252 d__2
= resid
, d__3
= (d__1
= b
[i__
- 1], abs(d__1
));
253 resid
= (max
)(d__2
, d__3
);
255 d__2
= normx
, d__3
= (d__1
= x
[i__
- 1], abs(d__1
));
256 normx
= (max
)(d__2
, d__3
);
259 eps
= epslon_(&c_b7
);
260 residn
= resid
/ (n
* norma
* normx
* eps
);
265 do_fio(&c__1, (char *)&residn, (ftnlen)sizeof(real_type));
266 do_fio(&c__1, (char *)&resid, (ftnlen)sizeof(real_type));
267 do_fio(&c__1, (char *)&eps, (ftnlen)sizeof(real_type));
268 do_fio(&c__1, (char *)&x[0], (ftnlen)sizeof(real_type));
269 do_fio(&c__1, (char *)&x[n - 1], (ftnlen)sizeof(real_type));
271 std::cout
<< std::setw(12) << std::setprecision(5) << residn
<< " " << resid
<< " " << eps
<< " " << x
[0] << " " << x
[n
- 1] << std::endl
;
275 do_fio(&c__1
, (char*)&n
, (ftnlen
)sizeof(integer
));
281 time
[3] = ops
/ (total
* 1e6
);
282 time
[4] = 2. / time
[3];
283 time
[5] = total
/ cray
;
285 do_fio(&c__1
, (char*)&lda
, (ftnlen
)sizeof(integer
));
288 for (i__
= 1; i__
<= 6; ++i__
)
290 // do_fio(&c__1, (char *)&time[i__ - 1], (ftnlen)sizeof(real_type));
291 std::cout
<< std::setw(12) << std::setprecision(5) << time
[i__
- 1];
295 do_lio(&c__9
, &c__1
, " end of tests -- this version dated 10/12/92", (ftnlen
)44);
298 s_stop("", (ftnlen
)0);
302 /* Subroutine */ int matgen_(real_type
* a
, integer
* lda
, integer
* n
,
303 real_type
* b
, real_type
* norma
)
305 /* System generated locals */
306 integer a_dim1
, a_offset
, i__1
, i__2
;
307 real_type d__1
, d__2
;
309 /* Local variables */
310 static integer i__
, j
;
311 static integer init
[4];
313 /* Parameter adjustments */
315 a_offset
= 1 + a_dim1
;
324 *norma
= CAST_TO_RT(0);
326 for (j
= 1; j
<= i__1
; ++j
)
329 for (i__
= 1; i__
<= i__2
; ++i__
)
331 a
[i__
+ j
* a_dim1
] = ran_(init
) - .5f
;
333 d__2
= (d__1
= a
[i__
+ j
* a_dim1
], abs(d__1
));
334 *norma
= (max
)(d__2
, *norma
);
340 for (i__
= 1; i__
<= i__1
; ++i__
)
342 b
[i__
] = CAST_TO_RT(0);
346 for (j
= 1; j
<= i__1
; ++j
)
349 for (i__
= 1; i__
<= i__2
; ++i__
)
351 b
[i__
] += a
[i__
+ j
* a_dim1
];
359 /* Subroutine */ int dgefa_(real_type
* a
, integer
* lda
, integer
* n
, integer
* ipvt
, integer
* info
)
361 /* System generated locals */
362 integer a_dim1
, a_offset
, i__1
, i__2
, i__3
;
364 /* Local variables */
365 static integer j
, k
, l
;
367 static integer kp1
, nm1
;
369 /* dgefa factors a double precision matrix by gaussian elimination. */
371 /* dgefa is usually called by dgeco, but it can be called */
372 /* directly with a saving in time if rcond is not needed. */
373 /* (time for dgeco) = (1 + 9/n)*(time for dgefa) . */
377 /* a double precision(lda, n) */
378 /* the matrix to be factored. */
381 /* the leading dimension of the array a . */
384 /* the order of the matrix a . */
388 /* a an upper triangular matrix and the multipliers */
389 /* which were used to obtain it. */
390 /* the factorization can be written a = l*u where */
391 /* l is a product of permutation and unit lower */
392 /* triangular matrices and u is upper triangular. */
394 /* ipvt integer(n) */
395 /* an integer vector of pivot indices. */
398 /* = 0 normal value. */
399 /* = k if u(k,k) .eq. 0.0 . this is not an error */
400 /* condition for this subroutine, but it does */
401 /* indicate that dgesl or dgedi will divide by zero */
402 /* if called. use rcond in dgeco for a reliable */
403 /* indication of singularity. */
405 /* linpack. this version dated 08/14/78 . */
406 /* cleve moler, university of new mexico, argonne national lab. */
408 /* subroutines and functions */
410 /* blas daxpy,dscal,idamax */
412 /* internal variables */
414 /* gaussian elimination with partial pivoting */
416 /* Parameter adjustments */
418 a_offset
= 1 + a_dim1
;
430 for (k
= 1; k
<= i__1
; ++k
)
434 /* find l = pivot index */
437 l
= idamax_(&i__2
, &a
[k
+ k
* a_dim1
], &c__1
) + k
- 1;
440 /* zero pivot implies this column already triangularized */
442 if (a
[l
+ k
* a_dim1
] == 0.)
447 /* interchange if necessary */
453 t
= a
[l
+ k
* a_dim1
];
454 a
[l
+ k
* a_dim1
] = a
[k
+ k
* a_dim1
];
455 a
[k
+ k
* a_dim1
] = t
;
458 /* compute multipliers */
460 t
= -1. / a
[k
+ k
* a_dim1
];
462 dscal_(&i__2
, &t
, &a
[k
+ 1 + k
* a_dim1
], &c__1
);
464 /* row elimination with column indexing */
467 for (j
= kp1
; j
<= i__2
; ++j
)
469 t
= a
[l
+ j
* a_dim1
];
474 a
[l
+ j
* a_dim1
] = a
[k
+ j
* a_dim1
];
475 a
[k
+ j
* a_dim1
] = t
;
478 daxpy_(&i__3
, &t
, &a
[k
+ 1 + k
* a_dim1
], &c__1
, &a
[k
+ 1 + j
* a_dim1
], &c__1
);
490 if (a
[*n
+ *n
* a_dim1
] == 0.)
497 /* Subroutine */ int dgesl_(real_type
* a
, integer
* lda
, integer
* n
, integer
* ipvt
, real_type
* b
, integer
* job
)
499 /* System generated locals */
500 integer a_dim1
, a_offset
, i__1
, i__2
;
502 /* Local variables */
505 static integer kb
, nm1
;
507 /* dgesl solves the double precision system */
508 /* a * x = b or trans(a) * x = b */
509 /* using the factors computed by dgeco or dgefa. */
513 /* a double precision(lda, n) */
514 /* the output from dgeco or dgefa. */
517 /* the leading dimension of the array a . */
520 /* the order of the matrix a . */
522 /* ipvt integer(n) */
523 /* the pivot vector from dgeco or dgefa. */
525 /* b double precision(n) */
526 /* the right hand side vector. */
529 /* = 0 to solve a*x = b , */
530 /* = nonzero to solve trans(a)*x = b where */
531 /* trans(a) is the transpose. */
535 /* b the solution vector x . */
537 /* error condition */
539 /* a division by zero will occur if the input factor contains a */
540 /* zero on the diagonal. technically this indicates singularity */
541 /* but it is often caused by improper arguments or improper */
542 /* setting of lda . it will not occur if the subroutines are */
543 /* called correctly and if dgeco has set rcond .gt. 0.0 */
544 /* or dgefa has set info .eq. 0 . */
546 /* to compute inverse(a) * c where c is a matrix */
548 /* call dgeco(a,lda,n,ipvt,rcond,z) */
549 /* if (rcond is too small) go to ... */
551 /* call dgesl(a,lda,n,ipvt,c(1,j),0) */
554 /* linpack. this version dated 08/14/78 . */
555 /* cleve moler, university of new mexico, argonne national lab. */
557 /* subroutines and functions */
559 /* blas daxpy,ddot */
561 /* internal variables */
563 /* Parameter adjustments */
565 a_offset
= 1 + a_dim1
;
577 /* job = 0 , solve a * x = b */
578 /* first solve l*y = b */
585 for (k
= 1; k
<= i__1
; ++k
)
597 daxpy_(&i__2
, &t
, &a
[k
+ 1 + k
* a_dim1
], &c__1
, &b
[k
+ 1], &c__1
);
602 /* now solve u*x = y */
605 for (kb
= 1; kb
<= i__1
; ++kb
)
608 b
[k
] /= a
[k
+ k
* a_dim1
];
611 daxpy_(&i__2
, &t
, &a
[k
* a_dim1
+ 1], &c__1
, &b
[1], &c__1
);
617 /* job = nonzero, solve trans(a) * x = b */
618 /* first solve trans(u)*y = b */
621 for (k
= 1; k
<= i__1
; ++k
)
624 t
= ddot_(&i__2
, &a
[k
* a_dim1
+ 1], &c__1
, &b
[1], &c__1
);
625 b
[k
] = (b
[k
] - t
) / a
[k
+ k
* a_dim1
];
629 /* now solve trans(l)*x = y */
636 for (kb
= 1; kb
<= i__1
; ++kb
)
640 b
[k
] += ddot_(&i__2
, &a
[k
+ 1 + k
* a_dim1
], &c__1
, &b
[k
+ 1], &c__1
);
658 /* Subroutine */ int daxpy_(integer
* n
, real_type
* da
, real_type
* dx
,
659 integer
* incx
, real_type
* dy
, integer
* incy
)
661 /* System generated locals */
664 /* Local variables */
665 static integer i__
, m
, ix
, iy
, mp1
;
667 /* constant times a vector plus a vector. */
668 /* uses unrolled loops for increments equal to one. */
669 /* jack dongarra, linpack, 3/11/78. */
671 /* Parameter adjustments */
684 if (*incx
== 1 && *incy
== 1)
689 /* code for unequal increments or equal increments */
696 ix
= (-(*n
) + 1) * *incx
+ 1;
700 iy
= (-(*n
) + 1) * *incy
+ 1;
703 for (i__
= 1; i__
<= i__1
; ++i__
)
705 dy
[iy
] += *da
* dx
[ix
];
712 /* code for both increments equal to 1 */
723 for (i__
= 1; i__
<= i__1
; ++i__
)
725 dy
[i__
] += *da
* dx
[i__
];
735 for (i__
= mp1
; i__
<= i__1
; i__
+= 4)
737 dy
[i__
] += *da
* dx
[i__
];
738 dy
[i__
+ 1] += *da
* dx
[i__
+ 1];
739 dy
[i__
+ 2] += *da
* dx
[i__
+ 2];
740 dy
[i__
+ 3] += *da
* dx
[i__
+ 3];
746 real_type
ddot_(integer
* n
, real_type
* dx
, integer
* incx
, real_type
* dy
,
749 /* System generated locals */
753 /* Local variables */
754 static integer i__
, m
, ix
, iy
, mp1
;
755 static real_type dtemp
;
757 /* forms the dot product of two vectors. */
758 /* uses unrolled loops for increments equal to one. */
759 /* jack dongarra, linpack, 3/11/78. */
761 /* Parameter adjustments */
766 ret_val
= CAST_TO_RT(0);
767 dtemp
= CAST_TO_RT(0);
772 if (*incx
== 1 && *incy
== 1)
777 /* code for unequal increments or equal increments */
784 ix
= (-(*n
) + 1) * *incx
+ 1;
788 iy
= (-(*n
) + 1) * *incy
+ 1;
791 for (i__
= 1; i__
<= i__1
; ++i__
)
793 dtemp
+= dx
[ix
] * dy
[iy
];
801 /* code for both increments equal to 1 */
812 for (i__
= 1; i__
<= i__1
; ++i__
)
814 dtemp
+= dx
[i__
] * dy
[i__
];
824 for (i__
= mp1
; i__
<= i__1
; i__
+= 5)
826 dtemp
= dtemp
+ dx
[i__
] * dy
[i__
] + dx
[i__
+ 1] * dy
[i__
+ 1] + dx
[i__
+ 2] * dy
[i__
+ 2] + dx
[i__
+ 3] * dy
[i__
+ 3] + dx
[i__
+ 4] * dy
[i__
+ 4];
834 /* Subroutine */ int dscal_(integer
* n
, real_type
* da
, real_type
* dx
,
837 /* System generated locals */
840 /* Local variables */
841 static integer i__
, m
, mp1
, nincx
;
843 /* scales a vector by a constant. */
844 /* uses unrolled loops for increment equal to one. */
845 /* jack dongarra, linpack, 3/11/78. */
847 /* Parameter adjustments */
860 /* code for increment not equal to 1 */
865 for (i__
= 1; i__2
< 0 ? i__
>= i__1
: i__
<= i__1
; i__
+= i__2
)
867 dx
[i__
] = *da
* dx
[i__
];
872 /* code for increment equal to 1 */
883 for (i__
= 1; i__
<= i__2
; ++i__
)
885 dx
[i__
] = *da
* dx
[i__
];
895 for (i__
= mp1
; i__
<= i__2
; i__
+= 5)
897 dx
[i__
] = *da
* dx
[i__
];
898 dx
[i__
+ 1] = *da
* dx
[i__
+ 1];
899 dx
[i__
+ 2] = *da
* dx
[i__
+ 2];
900 dx
[i__
+ 3] = *da
* dx
[i__
+ 3];
901 dx
[i__
+ 4] = *da
* dx
[i__
+ 4];
907 integer
idamax_(integer
* n
, real_type
* dx
, integer
* incx
)
909 /* System generated locals */
910 integer ret_val
, i__1
;
913 /* Local variables */
914 static integer i__
, ix
;
915 static real_type dmax__
;
917 /* finds the index of element having max. dabsolute value. */
918 /* jack dongarra, linpack, 3/11/78. */
920 /* Parameter adjustments */
939 /* code for increment not equal to 1 */
945 for (i__
= 2; i__
<= i__1
; ++i__
)
947 if ((d__1
= dx
[ix
], abs(d__1
)) <= dmax__
)
952 dmax__
= (d__1
= dx
[ix
], abs(d__1
));
959 /* code for increment equal to 1 */
964 for (i__
= 2; i__
<= i__1
; ++i__
)
966 if ((d__1
= dx
[i__
], abs(d__1
)) <= dmax__
)
971 dmax__
= (d__1
= dx
[i__
], abs(d__1
));
977 real_type
epslon_(real_type
* x
)
979 #if defined(TEST_MPF_100) || defined(TEST_MPFR_100) || defined(TEST_GMPXX) || defined(TEST_MPFRXX)
980 return std::ldexp(1.0, 1 - ((100 + 1) * 1000L) / 301L);
981 #elif defined(TEST_CPP_DEC_FLOAT_BN)
982 return std::pow(10.0, 1 - std::numeric_limits
<efx::cpp_dec_float_50
>::digits10
);
983 #elif defined(NATIVE_FLOAT128)
984 return FLT128_EPSILON
;
986 return CAST_TO_RT(std::numeric_limits
<real_type
>::epsilon());
990 /* Subroutine */ int mm_(real_type
* a
, integer
* lda
, integer
* n1
, integer
* n3
, real_type
* b
, integer
* ldb
, integer
* n2
, real_type
* c__
,
993 /* System generated locals */
994 integer a_dim1
, a_offset
, b_dim1
, b_offset
, c_dim1
, c_offset
, i__1
, i__2
;
996 /* Local variables */
997 static integer i__
, j
;
1000 /* multiply matrix b times matrix c and store the result in matrix a. */
1004 /* a double precision(lda,n3), matrix of n1 rows and n3 columns */
1006 /* lda integer, leading dimension of array a */
1008 /* n1 integer, number of rows in matrices a and b */
1010 /* n3 integer, number of columns in matrices a and c */
1012 /* b double precision(ldb,n2), matrix of n1 rows and n2 columns */
1014 /* ldb integer, leading dimension of array b */
1016 /* n2 integer, number of columns in matrix b, and number of rows in */
1019 /* c double precision(ldc,n3), matrix of n2 rows and n3 columns */
1021 /* ldc integer, leading dimension of array c */
1023 /* ---------------------------------------------------------------------- */
1025 /* Parameter adjustments */
1027 a_offset
= 1 + a_dim1
;
1030 b_offset
= 1 + b_dim1
;
1033 c_offset
= 1 + c_dim1
;
1038 for (j
= 1; j
<= i__1
; ++j
)
1041 for (i__
= 1; i__
<= i__2
; ++i__
)
1043 a
[i__
+ j
* a_dim1
] = CAST_TO_RT(0);
1046 dmxpy_(n2
, &a
[j
* a_dim1
+ 1], n1
, ldb
, &c__
[j
* c_dim1
+ 1], &b
[b_offset
]);
1053 /* Subroutine */ int dmxpy_(integer
* n1
, real_type
* y
, integer
* n2
, integer
* ldm
, real_type
* x
, real_type
* m
)
1055 /* System generated locals */
1056 integer m_dim1
, m_offset
, i__1
, i__2
;
1058 /* Local variables */
1059 static integer i__
, j
, jmin
;
1062 /* multiply matrix m times vector x and add the result to vector y. */
1066 /* n1 integer, number of elements in vector y, and number of rows in */
1069 /* y double precision(n1), vector of length n1 to which is added */
1070 /* the product m*x */
1072 /* n2 integer, number of elements in vector x, and number of columns */
1075 /* ldm integer, leading dimension of array m */
1077 /* x double precision(n2), vector of length n2 */
1079 /* m double precision(ldm,n2), matrix of n1 rows and n2 columns */
1081 /* ---------------------------------------------------------------------- */
1083 /* cleanup odd vector */
1085 /* Parameter adjustments */
1088 m_offset
= 1 + m_dim1
;
1097 for (i__
= 1; i__
<= i__1
; ++i__
)
1099 y
[i__
] += x
[j
] * m
[i__
+ j
* m_dim1
];
1104 /* cleanup odd group of two vectors */
1110 for (i__
= 1; i__
<= i__1
; ++i__
)
1112 y
[i__
] = y
[i__
] + x
[j
- 1] * m
[i__
+ (j
- 1) * m_dim1
] + x
[j
] * m
[i__
+ j
* m_dim1
];
1117 /* cleanup odd group of four vectors */
1123 for (i__
= 1; i__
<= i__1
; ++i__
)
1125 y
[i__
] = y
[i__
] + x
[j
- 3] * m
[i__
+ (j
- 3) * m_dim1
] + x
[j
- 2] * m
[i__
+ (j
- 2) * m_dim1
] + x
[j
- 1] * m
[i__
+ (j
- 1) * m_dim1
] + x
[j
] * m
[i__
+ j
* m_dim1
];
1130 /* cleanup odd group of eight vectors */
1136 for (i__
= 1; i__
<= i__1
; ++i__
)
1138 y
[i__
] = y
[i__
] + x
[j
- 7] * m
[i__
+ (j
- 7) * m_dim1
] + x
[j
- 6] * m
[i__
+ (j
- 6) * m_dim1
] + x
[j
- 5] * m
[i__
+ (j
- 5) * m_dim1
] + x
[j
- 4] * m
[i__
+ (j
- 4) * m_dim1
] + x
[j
- 3] * m
[i__
+ (j
- 3) * m_dim1
] + x
[j
- 2] * m
[i__
+ (j
- 2) * m_dim1
] + x
[j
- 1] * m
[i__
+ (j
- 1) * m_dim1
] + x
[j
] * m
[i__
+ j
* m_dim1
];
1143 /* main loop - groups of sixteen vectors */
1147 for (j
= jmin
; j
<= i__1
; j
+= 16)
1150 for (i__
= 1; i__
<= i__2
; ++i__
)
1152 y
[i__
] = y
[i__
] + x
[j
- 15] * m
[i__
+ (j
- 15) * m_dim1
] + x
[j
- 14] * m
[i__
+ (j
- 14) * m_dim1
] + x
[j
- 13] * m
[i__
+ (j
- 13) * m_dim1
] + x
[j
- 12] * m
[i__
+ (j
- 12) * m_dim1
] + x
[j
- 11] * m
[i__
+ (j
- 11) * m_dim1
] + x
[j
- 10] * m
[i__
+ (j
- 10) * m_dim1
] + x
[j
- 9] * m
[i__
+ (j
- 9) * m_dim1
] + x
[j
- 8] * m
[i__
+ (j
- 8) * m_dim1
] + x
[j
- 7] * m
[i__
+ (j
- 7) * m_dim1
] + x
[j
- 6] * m
[i__
+ (j
- 6) * m_dim1
] + x
[j
- 5] * m
[i__
+ (j
- 5) * m_dim1
] + x
[j
- 4] * m
[i__
+ (j
- 4) * m_dim1
] + x
[j
- 3] * m
[i__
+ (j
- 3) * m_dim1
] + x
[j
- 2] * m
[i__
+ (j
- 2) * m_dim1
] + x
[j
- 1] * m
[i__
+ (j
- 1) * m_dim1
] + x
[j
] * m
[i__
+ j
* m_dim1
];
1160 real_type
ran_(integer
* iseed
)
1162 /* System generated locals */
1165 /* Local variables */
1166 static integer it1
, it2
, it3
, it4
;
1168 /* modified from the LAPACK auxiliary routine 10/12/92 JD */
1169 /* -- LAPACK auxiliary routine (version 1.0) -- */
1170 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
1171 /* Courant Institute, Argonne National Lab, and Rice University */
1172 /* February 29, 1992 */
1174 /* .. Array Arguments .. */
1180 /* DLARAN returns a random double number from a uniform (0,1) */
1186 /* ISEED (input/output) INTEGER array, dimension (4) */
1187 /* On entry, the seed of the random number generator; the array */
1188 /* elements must be between 0 and 4095, and ISEED(4) must be */
1190 /* On exit, the seed is updated. */
1192 /* Further Details */
1193 /* =============== */
1195 /* This routine uses a multiplicative congruential method with modulus */
1196 /* 2**48 and multiplier 33952834046453 (see G.S.Fishman, */
1197 /* 'Multiplicative congruential random number generators with modulus */
1198 /* 2**b: an exhaustive analysis for b = 32 and a partial analysis for */
1199 /* b = 48', Math. Comp. 189, pp 331-344, 1990). */
1201 /* 48-bit integers are stored in 4 integer array elements with 12 bits */
1202 /* per element. Hence the routine is portable across machines with */
1203 /* integers of 32 bits or more. */
1205 /* .. Parameters .. */
1207 /* .. Local Scalars .. */
1209 /* .. Intrinsic Functions .. */
1211 /* .. Executable Statements .. */
1213 /* multiply the seed by the multiplier modulo 2**48 */
1215 /* Parameter adjustments */
1219 it4
= iseed
[4] * 2549;
1222 it3
= it3
+ iseed
[3] * 2549 + iseed
[4] * 2508;
1225 it2
= it2
+ iseed
[2] * 2549 + iseed
[3] * 2508 + iseed
[4] * 322;
1228 it1
= it1
+ iseed
[1] * 2549 + iseed
[2] * 2508 + iseed
[3] * 322 + iseed
[4] * 494;
1231 /* return updated seed */
1238 /* convert 48-bit integer to a double number in the interval (0,1) */
1240 ret_val
= ((real_type
)it1
+ ((real_type
)it2
+ ((real_type
)it3
+ (real_type
)it4
* 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4;
1252 norm. resid resid machep x(1) x(n)
1253 6.4915 7.207e-013 2.2204e-016 1 1
1257 times are reported for matrices of order 1000
1258 factor solve total mflops unit ratio
1259 times for array with leading dimension of1001
1260 1.443 0.003 1.446 462.43 0.004325 25.821
1266 norm. resid resid machep x(1) x(n)
1267 3.6575e-05 5.2257e-103 2.8575e-101 1 1
1271 times are reported for matrices of order 1000
1272 factor solve total mflops unit ratio
1273 times for array with leading dimension of1001
1274 266.45 0.798 267.24 2.5021 0.79933 4772.2
1277 number<gmp_float<100> >:
1278 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1280 norm. resid resid machep x(1) x(n)
1281 0.36575e-4 0.52257e-102 0.28575e-100 0.1e1 0.1e1
1285 times are reported for matrices of order 1000
1286 factor solve total mflops unit ratio
1287 times for array with leading dimension of1001
1288 279.96 0.84 280.8 2.3813 0.83988 5014.3
1290 boost::multiprecision::ef::cpp_dec_float_50:
1291 ~~~~~~~~~~~~~~~~~~~~~~~~~
1293 norm. resid resid machep x(1) x(n)
1294 2.551330735e-16 1.275665107e-112 1e-99 1 1
1298 times are reported for matrices of order 1000
1299 factor solve total mflops unit ratio
1300 times for array with leading dimension of1001
1301 363.89 1.074 364.97 1.8321 1.0916 6517.3