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)
5
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
12 cc *.o -lf2c -lm
13 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
14
15 http://www.netlib.org/f2c/libf2c.zip
16 */
17 #include <iostream>
18 #include <iomanip>
19 #include <cmath>
20
21 #if defined(TEST_GMPXX)
22 #include <gmpxx.h>
23 typedef mpf_class real_type;
24 #elif defined(TEST_MPFRXX)
25 #include <gmpfrxx.h>
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;
39
operator <<(std::ostream & os,const __float128 & f)40 std::ostream& operator<<(std::ostream& os, const __float128& f)
41 {
42 return os << boost::multiprecision::float128(f);
43 }
44
45 #include <boost/type_traits/has_left_shift.hpp>
46
47 namespace boost {
48
49 template <>
50 struct has_left_shift<std::basic_ostream<char>, __float128> : public mpl::true_
51 {};
52
53 template <>
lexical_cast(const __float128 & f)54 double lexical_cast<double, __float128>(const __float128& f)
55 {
56 return f;
57 }
58
59 } // namespace boost
60
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;
70 #else
71 typedef double real_type;
72 #endif
73
74 #include <boost/lexical_cast.hpp>
75
76 #ifndef CAST_TO_RT
77 #define CAST_TO_RT(x) x
78 #endif
79
80 extern "C" {
81 #include "f2c.h"
82 integer s_wsfe(cilist*), e_wsfe(void), do_fio(integer*, char*, ftnlen),
83 s_wsle(cilist*), do_lio(integer*, integer*, char*, ftnlen),
84 e_wsle(void);
85 /* Subroutine */ int s_stop(char*, ftnlen);
86
87 #undef abs
88 #undef dabs
89 #define dabs abs
90 #undef dmin
91 #undef dmax
92 #define dmin min
93 #define dmax max
94 }
95 #include <time.h>
96
97 using std::max;
98 using std::min;
99
100 /* Table of constant values */
101
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;
106
second_(void)107 inline double second_(void)
108 {
109 return ((double)(clock())) / CLOCKS_PER_SEC;
110 }
111
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*);
123
MAIN__()124 extern "C" int MAIN__()
125 {
126 #ifdef TEST_MPF_50
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;
142 #else
143 std::cout << "Testing double" << std::endl;
144 #endif
145
146 /* Format strings */
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"
151 "s.utk.edu\002/)";
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"
156 "der \002,i5)";
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"
160 "f\002,i4)";
161 static char fmt_110[] = "(6(1pe11.3))";
162
163 /* System generated locals */
164 integer i__1;
165 real_type d__1, d__2, d__3;
166
167 /* Builtin functions */
168
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];
173 static double t1;
174 static integer lda;
175 static double ops;
176 static real_type eps;
177 static integer info;
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;
183
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};
193
194 lda = 1001;
195
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. */
200
201 n = 1000;
202 cray = .056f;
203 s_wsfe(&io___4);
204 e_wsfe();
205 /* Computing 3rd power */
206 d__1 = (real_type)n;
207 /* Computing 2nd power */
208 d__2 = (real_type)n;
209 ops = boost::lexical_cast<double>(real_type(d__1 * (d__1 * d__1) * 2. / 3. + d__2 * d__2 * 2.));
210
211 matgen_(a, &lda, &n, b, &norma);
212
213 /* ****************************************************************** */
214 /* ****************************************************************** */
215 /* you should replace the call to dgefa and dgesl */
216 /* by calls to your linear equation solver. */
217 /* ****************************************************************** */
218 /* ****************************************************************** */
219
220 t1 = second_();
221 dgefa_(a, &lda, &n, ipvt, &info);
222 time[0] = second_() - t1;
223 t1 = second_();
224 dgesl_(a, &lda, &n, ipvt, b, &c__0);
225 time[1] = second_() - t1;
226 total = time[0] + time[1];
227 /* ****************************************************************** */
228 /* ****************************************************************** */
229
230 /* compute a residual to verify results. */
231
232 i__1 = n;
233 for (i__ = 1; i__ <= i__1; ++i__)
234 {
235 x[i__ - 1] = b[i__ - 1];
236 /* L10: */
237 }
238 matgen_(a, &lda, &n, b, &norma);
239 i__1 = n;
240 for (i__ = 1; i__ <= i__1; ++i__)
241 {
242 b[i__ - 1] = -b[i__ - 1];
243 /* L20: */
244 }
245 dmxpy_(&n, b, &n, &lda, x, a);
246 resid = CAST_TO_RT(0);
247 normx = CAST_TO_RT(0);
248 i__1 = n;
249 for (i__ = 1; i__ <= i__1; ++i__)
250 {
251 /* Computing MAX */
252 d__2 = resid, d__3 = (d__1 = b[i__ - 1], abs(d__1));
253 resid = (max)(d__2, d__3);
254 /* Computing MAX */
255 d__2 = normx, d__3 = (d__1 = x[i__ - 1], abs(d__1));
256 normx = (max)(d__2, d__3);
257 /* L30: */
258 }
259 eps = epslon_(&c_b7);
260 residn = resid / (n * norma * normx * eps);
261 s_wsfe(&io___20);
262 e_wsfe();
263 s_wsfe(&io___21);
264 /*
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));
270 */
271 std::cout << std::setw(12) << std::setprecision(5) << residn << " " << resid << " " << eps << " " << x[0] << " " << x[n - 1] << std::endl;
272 e_wsfe();
273
274 s_wsfe(&io___22);
275 do_fio(&c__1, (char*)&n, (ftnlen)sizeof(integer));
276 e_wsfe();
277 s_wsfe(&io___23);
278 e_wsfe();
279
280 time[2] = total;
281 time[3] = ops / (total * 1e6);
282 time[4] = 2. / time[3];
283 time[5] = total / cray;
284 s_wsfe(&io___24);
285 do_fio(&c__1, (char*)&lda, (ftnlen)sizeof(integer));
286 e_wsfe();
287 s_wsfe(&io___25);
288 for (i__ = 1; i__ <= 6; ++i__)
289 {
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];
292 }
293 e_wsfe();
294 s_wsle(&io___26);
295 do_lio(&c__9, &c__1, " end of tests -- this version dated 10/12/92", (ftnlen)44);
296 e_wsle();
297
298 s_stop("", (ftnlen)0);
299 return 0;
300 } /* MAIN__ */
301
matgen_(real_type * a,integer * lda,integer * n,real_type * b,real_type * norma)302 /* Subroutine */ int matgen_(real_type* a, integer* lda, integer* n,
303 real_type* b, real_type* norma)
304 {
305 /* System generated locals */
306 integer a_dim1, a_offset, i__1, i__2;
307 real_type d__1, d__2;
308
309 /* Local variables */
310 static integer i__, j;
311 static integer init[4];
312
313 /* Parameter adjustments */
314 a_dim1 = *lda;
315 a_offset = 1 + a_dim1;
316 a -= a_offset;
317 --b;
318
319 /* Function Body */
320 init[0] = 1;
321 init[1] = 2;
322 init[2] = 3;
323 init[3] = 1325;
324 *norma = CAST_TO_RT(0);
325 i__1 = *n;
326 for (j = 1; j <= i__1; ++j)
327 {
328 i__2 = *n;
329 for (i__ = 1; i__ <= i__2; ++i__)
330 {
331 a[i__ + j * a_dim1] = ran_(init) - .5f;
332 /* Computing MAX */
333 d__2 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
334 *norma = (max)(d__2, *norma);
335 /* L20: */
336 }
337 /* L30: */
338 }
339 i__1 = *n;
340 for (i__ = 1; i__ <= i__1; ++i__)
341 {
342 b[i__] = CAST_TO_RT(0);
343 /* L35: */
344 }
345 i__1 = *n;
346 for (j = 1; j <= i__1; ++j)
347 {
348 i__2 = *n;
349 for (i__ = 1; i__ <= i__2; ++i__)
350 {
351 b[i__] += a[i__ + j * a_dim1];
352 /* L40: */
353 }
354 /* L50: */
355 }
356 return 0;
357 } /* matgen_ */
358
dgefa_(real_type * a,integer * lda,integer * n,integer * ipvt,integer * info)359 /* Subroutine */ int dgefa_(real_type* a, integer* lda, integer* n, integer* ipvt, integer* info)
360 {
361 /* System generated locals */
362 integer a_dim1, a_offset, i__1, i__2, i__3;
363
364 /* Local variables */
365 static integer j, k, l;
366 static real_type t;
367 static integer kp1, nm1;
368
369 /* dgefa factors a double precision matrix by gaussian elimination. */
370
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) . */
374
375 /* on entry */
376
377 /* a double precision(lda, n) */
378 /* the matrix to be factored. */
379
380 /* lda integer */
381 /* the leading dimension of the array a . */
382
383 /* n integer */
384 /* the order of the matrix a . */
385
386 /* on return */
387
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. */
393
394 /* ipvt integer(n) */
395 /* an integer vector of pivot indices. */
396
397 /* info integer */
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. */
404
405 /* linpack. this version dated 08/14/78 . */
406 /* cleve moler, university of new mexico, argonne national lab. */
407
408 /* subroutines and functions */
409
410 /* blas daxpy,dscal,idamax */
411
412 /* internal variables */
413
414 /* gaussian elimination with partial pivoting */
415
416 /* Parameter adjustments */
417 a_dim1 = *lda;
418 a_offset = 1 + a_dim1;
419 a -= a_offset;
420 --ipvt;
421
422 /* Function Body */
423 *info = 0;
424 nm1 = *n - 1;
425 if (nm1 < 1)
426 {
427 goto L70;
428 }
429 i__1 = nm1;
430 for (k = 1; k <= i__1; ++k)
431 {
432 kp1 = k + 1;
433
434 /* find l = pivot index */
435
436 i__2 = *n - k + 1;
437 l = idamax_(&i__2, &a[k + k * a_dim1], &c__1) + k - 1;
438 ipvt[k] = l;
439
440 /* zero pivot implies this column already triangularized */
441
442 if (a[l + k * a_dim1] == 0.)
443 {
444 goto L40;
445 }
446
447 /* interchange if necessary */
448
449 if (l == k)
450 {
451 goto L10;
452 }
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;
456 L10:
457
458 /* compute multipliers */
459
460 t = -1. / a[k + k * a_dim1];
461 i__2 = *n - k;
462 dscal_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1);
463
464 /* row elimination with column indexing */
465
466 i__2 = *n;
467 for (j = kp1; j <= i__2; ++j)
468 {
469 t = a[l + j * a_dim1];
470 if (l == k)
471 {
472 goto L20;
473 }
474 a[l + j * a_dim1] = a[k + j * a_dim1];
475 a[k + j * a_dim1] = t;
476 L20:
477 i__3 = *n - k;
478 daxpy_(&i__3, &t, &a[k + 1 + k * a_dim1], &c__1, &a[k + 1 + j * a_dim1], &c__1);
479 /* L30: */
480 }
481 goto L50;
482 L40:
483 *info = k;
484 L50:
485 /* L60: */
486 ;
487 }
488 L70:
489 ipvt[*n] = *n;
490 if (a[*n + *n * a_dim1] == 0.)
491 {
492 *info = *n;
493 }
494 return 0;
495 } /* dgefa_ */
496
dgesl_(real_type * a,integer * lda,integer * n,integer * ipvt,real_type * b,integer * job)497 /* Subroutine */ int dgesl_(real_type* a, integer* lda, integer* n, integer* ipvt, real_type* b, integer* job)
498 {
499 /* System generated locals */
500 integer a_dim1, a_offset, i__1, i__2;
501
502 /* Local variables */
503 static integer k, l;
504 static real_type t;
505 static integer kb, nm1;
506
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. */
510
511 /* on entry */
512
513 /* a double precision(lda, n) */
514 /* the output from dgeco or dgefa. */
515
516 /* lda integer */
517 /* the leading dimension of the array a . */
518
519 /* n integer */
520 /* the order of the matrix a . */
521
522 /* ipvt integer(n) */
523 /* the pivot vector from dgeco or dgefa. */
524
525 /* b double precision(n) */
526 /* the right hand side vector. */
527
528 /* job integer */
529 /* = 0 to solve a*x = b , */
530 /* = nonzero to solve trans(a)*x = b where */
531 /* trans(a) is the transpose. */
532
533 /* on return */
534
535 /* b the solution vector x . */
536
537 /* error condition */
538
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 . */
545
546 /* to compute inverse(a) * c where c is a matrix */
547 /* with p columns */
548 /* call dgeco(a,lda,n,ipvt,rcond,z) */
549 /* if (rcond is too small) go to ... */
550 /* do 10 j = 1, p */
551 /* call dgesl(a,lda,n,ipvt,c(1,j),0) */
552 /* 10 continue */
553
554 /* linpack. this version dated 08/14/78 . */
555 /* cleve moler, university of new mexico, argonne national lab. */
556
557 /* subroutines and functions */
558
559 /* blas daxpy,ddot */
560
561 /* internal variables */
562
563 /* Parameter adjustments */
564 a_dim1 = *lda;
565 a_offset = 1 + a_dim1;
566 a -= a_offset;
567 --ipvt;
568 --b;
569
570 /* Function Body */
571 nm1 = *n - 1;
572 if (*job != 0)
573 {
574 goto L50;
575 }
576
577 /* job = 0 , solve a * x = b */
578 /* first solve l*y = b */
579
580 if (nm1 < 1)
581 {
582 goto L30;
583 }
584 i__1 = nm1;
585 for (k = 1; k <= i__1; ++k)
586 {
587 l = ipvt[k];
588 t = b[l];
589 if (l == k)
590 {
591 goto L10;
592 }
593 b[l] = b[k];
594 b[k] = t;
595 L10:
596 i__2 = *n - k;
597 daxpy_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
598 /* L20: */
599 }
600 L30:
601
602 /* now solve u*x = y */
603
604 i__1 = *n;
605 for (kb = 1; kb <= i__1; ++kb)
606 {
607 k = *n + 1 - kb;
608 b[k] /= a[k + k * a_dim1];
609 t = -b[k];
610 i__2 = k - 1;
611 daxpy_(&i__2, &t, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
612 /* L40: */
613 }
614 goto L100;
615 L50:
616
617 /* job = nonzero, solve trans(a) * x = b */
618 /* first solve trans(u)*y = b */
619
620 i__1 = *n;
621 for (k = 1; k <= i__1; ++k)
622 {
623 i__2 = k - 1;
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];
626 /* L60: */
627 }
628
629 /* now solve trans(l)*x = y */
630
631 if (nm1 < 1)
632 {
633 goto L90;
634 }
635 i__1 = nm1;
636 for (kb = 1; kb <= i__1; ++kb)
637 {
638 k = *n - kb;
639 i__2 = *n - k;
640 b[k] += ddot_(&i__2, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
641 l = ipvt[k];
642 if (l == k)
643 {
644 goto L70;
645 }
646 t = b[l];
647 b[l] = b[k];
648 b[k] = t;
649 L70:
650 /* L80: */
651 ;
652 }
653 L90:
654 L100:
655 return 0;
656 } /* dgesl_ */
657
daxpy_(integer * n,real_type * da,real_type * dx,integer * incx,real_type * dy,integer * incy)658 /* Subroutine */ int daxpy_(integer* n, real_type* da, real_type* dx,
659 integer* incx, real_type* dy, integer* incy)
660 {
661 /* System generated locals */
662 integer i__1;
663
664 /* Local variables */
665 static integer i__, m, ix, iy, mp1;
666
667 /* constant times a vector plus a vector. */
668 /* uses unrolled loops for increments equal to one. */
669 /* jack dongarra, linpack, 3/11/78. */
670
671 /* Parameter adjustments */
672 --dy;
673 --dx;
674
675 /* Function Body */
676 if (*n <= 0)
677 {
678 return 0;
679 }
680 if (*da == 0.)
681 {
682 return 0;
683 }
684 if (*incx == 1 && *incy == 1)
685 {
686 goto L20;
687 }
688
689 /* code for unequal increments or equal increments */
690 /* not equal to 1 */
691
692 ix = 1;
693 iy = 1;
694 if (*incx < 0)
695 {
696 ix = (-(*n) + 1) * *incx + 1;
697 }
698 if (*incy < 0)
699 {
700 iy = (-(*n) + 1) * *incy + 1;
701 }
702 i__1 = *n;
703 for (i__ = 1; i__ <= i__1; ++i__)
704 {
705 dy[iy] += *da * dx[ix];
706 ix += *incx;
707 iy += *incy;
708 /* L10: */
709 }
710 return 0;
711
712 /* code for both increments equal to 1 */
713
714 /* clean-up loop */
715
716 L20:
717 m = *n % 4;
718 if (m == 0)
719 {
720 goto L40;
721 }
722 i__1 = m;
723 for (i__ = 1; i__ <= i__1; ++i__)
724 {
725 dy[i__] += *da * dx[i__];
726 /* L30: */
727 }
728 if (*n < 4)
729 {
730 return 0;
731 }
732 L40:
733 mp1 = m + 1;
734 i__1 = *n;
735 for (i__ = mp1; i__ <= i__1; i__ += 4)
736 {
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];
741 /* L50: */
742 }
743 return 0;
744 } /* daxpy_ */
745
ddot_(integer * n,real_type * dx,integer * incx,real_type * dy,integer * incy)746 real_type ddot_(integer* n, real_type* dx, integer* incx, real_type* dy,
747 integer* incy)
748 {
749 /* System generated locals */
750 integer i__1;
751 real_type ret_val;
752
753 /* Local variables */
754 static integer i__, m, ix, iy, mp1;
755 static real_type dtemp;
756
757 /* forms the dot product of two vectors. */
758 /* uses unrolled loops for increments equal to one. */
759 /* jack dongarra, linpack, 3/11/78. */
760
761 /* Parameter adjustments */
762 --dy;
763 --dx;
764
765 /* Function Body */
766 ret_val = CAST_TO_RT(0);
767 dtemp = CAST_TO_RT(0);
768 if (*n <= 0)
769 {
770 return ret_val;
771 }
772 if (*incx == 1 && *incy == 1)
773 {
774 goto L20;
775 }
776
777 /* code for unequal increments or equal increments */
778 /* not equal to 1 */
779
780 ix = 1;
781 iy = 1;
782 if (*incx < 0)
783 {
784 ix = (-(*n) + 1) * *incx + 1;
785 }
786 if (*incy < 0)
787 {
788 iy = (-(*n) + 1) * *incy + 1;
789 }
790 i__1 = *n;
791 for (i__ = 1; i__ <= i__1; ++i__)
792 {
793 dtemp += dx[ix] * dy[iy];
794 ix += *incx;
795 iy += *incy;
796 /* L10: */
797 }
798 ret_val = dtemp;
799 return ret_val;
800
801 /* code for both increments equal to 1 */
802
803 /* clean-up loop */
804
805 L20:
806 m = *n % 5;
807 if (m == 0)
808 {
809 goto L40;
810 }
811 i__1 = m;
812 for (i__ = 1; i__ <= i__1; ++i__)
813 {
814 dtemp += dx[i__] * dy[i__];
815 /* L30: */
816 }
817 if (*n < 5)
818 {
819 goto L60;
820 }
821 L40:
822 mp1 = m + 1;
823 i__1 = *n;
824 for (i__ = mp1; i__ <= i__1; i__ += 5)
825 {
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];
827 /* L50: */
828 }
829 L60:
830 ret_val = dtemp;
831 return ret_val;
832 } /* ddot_ */
833
dscal_(integer * n,real_type * da,real_type * dx,integer * incx)834 /* Subroutine */ int dscal_(integer* n, real_type* da, real_type* dx,
835 integer* incx)
836 {
837 /* System generated locals */
838 integer i__1, i__2;
839
840 /* Local variables */
841 static integer i__, m, mp1, nincx;
842
843 /* scales a vector by a constant. */
844 /* uses unrolled loops for increment equal to one. */
845 /* jack dongarra, linpack, 3/11/78. */
846
847 /* Parameter adjustments */
848 --dx;
849
850 /* Function Body */
851 if (*n <= 0)
852 {
853 return 0;
854 }
855 if (*incx == 1)
856 {
857 goto L20;
858 }
859
860 /* code for increment not equal to 1 */
861
862 nincx = *n * *incx;
863 i__1 = nincx;
864 i__2 = *incx;
865 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2)
866 {
867 dx[i__] = *da * dx[i__];
868 /* L10: */
869 }
870 return 0;
871
872 /* code for increment equal to 1 */
873
874 /* clean-up loop */
875
876 L20:
877 m = *n % 5;
878 if (m == 0)
879 {
880 goto L40;
881 }
882 i__2 = m;
883 for (i__ = 1; i__ <= i__2; ++i__)
884 {
885 dx[i__] = *da * dx[i__];
886 /* L30: */
887 }
888 if (*n < 5)
889 {
890 return 0;
891 }
892 L40:
893 mp1 = m + 1;
894 i__2 = *n;
895 for (i__ = mp1; i__ <= i__2; i__ += 5)
896 {
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];
902 /* L50: */
903 }
904 return 0;
905 } /* dscal_ */
906
idamax_(integer * n,real_type * dx,integer * incx)907 integer idamax_(integer* n, real_type* dx, integer* incx)
908 {
909 /* System generated locals */
910 integer ret_val, i__1;
911 real_type d__1;
912
913 /* Local variables */
914 static integer i__, ix;
915 static real_type dmax__;
916
917 /* finds the index of element having max. dabsolute value. */
918 /* jack dongarra, linpack, 3/11/78. */
919
920 /* Parameter adjustments */
921 --dx;
922
923 /* Function Body */
924 ret_val = 0;
925 if (*n < 1)
926 {
927 return ret_val;
928 }
929 ret_val = 1;
930 if (*n == 1)
931 {
932 return ret_val;
933 }
934 if (*incx == 1)
935 {
936 goto L20;
937 }
938
939 /* code for increment not equal to 1 */
940
941 ix = 1;
942 dmax__ = abs(dx[1]);
943 ix += *incx;
944 i__1 = *n;
945 for (i__ = 2; i__ <= i__1; ++i__)
946 {
947 if ((d__1 = dx[ix], abs(d__1)) <= dmax__)
948 {
949 goto L5;
950 }
951 ret_val = i__;
952 dmax__ = (d__1 = dx[ix], abs(d__1));
953 L5:
954 ix += *incx;
955 /* L10: */
956 }
957 return ret_val;
958
959 /* code for increment equal to 1 */
960
961 L20:
962 dmax__ = abs(dx[1]);
963 i__1 = *n;
964 for (i__ = 2; i__ <= i__1; ++i__)
965 {
966 if ((d__1 = dx[i__], abs(d__1)) <= dmax__)
967 {
968 goto L30;
969 }
970 ret_val = i__;
971 dmax__ = (d__1 = dx[i__], abs(d__1));
972 L30:;
973 }
974 return ret_val;
975 } /* idamax_ */
976
epslon_(real_type * x)977 real_type epslon_(real_type* x)
978 {
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;
985 #else
986 return CAST_TO_RT(std::numeric_limits<real_type>::epsilon());
987 #endif
988 } /* epslon_ */
989
mm_(real_type * a,integer * lda,integer * n1,integer * n3,real_type * b,integer * ldb,integer * n2,real_type * c__,integer * ldc)990 /* Subroutine */ int mm_(real_type* a, integer* lda, integer* n1, integer* n3, real_type* b, integer* ldb, integer* n2, real_type* c__,
991 integer* ldc)
992 {
993 /* System generated locals */
994 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2;
995
996 /* Local variables */
997 static integer i__, j;
998
999 /* purpose: */
1000 /* multiply matrix b times matrix c and store the result in matrix a. */
1001
1002 /* parameters: */
1003
1004 /* a double precision(lda,n3), matrix of n1 rows and n3 columns */
1005
1006 /* lda integer, leading dimension of array a */
1007
1008 /* n1 integer, number of rows in matrices a and b */
1009
1010 /* n3 integer, number of columns in matrices a and c */
1011
1012 /* b double precision(ldb,n2), matrix of n1 rows and n2 columns */
1013
1014 /* ldb integer, leading dimension of array b */
1015
1016 /* n2 integer, number of columns in matrix b, and number of rows in */
1017 /* matrix c */
1018
1019 /* c double precision(ldc,n3), matrix of n2 rows and n3 columns */
1020
1021 /* ldc integer, leading dimension of array c */
1022
1023 /* ---------------------------------------------------------------------- */
1024
1025 /* Parameter adjustments */
1026 a_dim1 = *lda;
1027 a_offset = 1 + a_dim1;
1028 a -= a_offset;
1029 b_dim1 = *ldb;
1030 b_offset = 1 + b_dim1;
1031 b -= b_offset;
1032 c_dim1 = *ldc;
1033 c_offset = 1 + c_dim1;
1034 c__ -= c_offset;
1035
1036 /* Function Body */
1037 i__1 = *n3;
1038 for (j = 1; j <= i__1; ++j)
1039 {
1040 i__2 = *n1;
1041 for (i__ = 1; i__ <= i__2; ++i__)
1042 {
1043 a[i__ + j * a_dim1] = CAST_TO_RT(0);
1044 /* L10: */
1045 }
1046 dmxpy_(n2, &a[j * a_dim1 + 1], n1, ldb, &c__[j * c_dim1 + 1], &b[b_offset]);
1047 /* L20: */
1048 }
1049
1050 return 0;
1051 } /* mm_ */
1052
dmxpy_(integer * n1,real_type * y,integer * n2,integer * ldm,real_type * x,real_type * m)1053 /* Subroutine */ int dmxpy_(integer* n1, real_type* y, integer* n2, integer* ldm, real_type* x, real_type* m)
1054 {
1055 /* System generated locals */
1056 integer m_dim1, m_offset, i__1, i__2;
1057
1058 /* Local variables */
1059 static integer i__, j, jmin;
1060
1061 /* purpose: */
1062 /* multiply matrix m times vector x and add the result to vector y. */
1063
1064 /* parameters: */
1065
1066 /* n1 integer, number of elements in vector y, and number of rows in */
1067 /* matrix m */
1068
1069 /* y double precision(n1), vector of length n1 to which is added */
1070 /* the product m*x */
1071
1072 /* n2 integer, number of elements in vector x, and number of columns */
1073 /* in matrix m */
1074
1075 /* ldm integer, leading dimension of array m */
1076
1077 /* x double precision(n2), vector of length n2 */
1078
1079 /* m double precision(ldm,n2), matrix of n1 rows and n2 columns */
1080
1081 /* ---------------------------------------------------------------------- */
1082
1083 /* cleanup odd vector */
1084
1085 /* Parameter adjustments */
1086 --y;
1087 m_dim1 = *ldm;
1088 m_offset = 1 + m_dim1;
1089 m -= m_offset;
1090 --x;
1091
1092 /* Function Body */
1093 j = *n2 % 2;
1094 if (j >= 1)
1095 {
1096 i__1 = *n1;
1097 for (i__ = 1; i__ <= i__1; ++i__)
1098 {
1099 y[i__] += x[j] * m[i__ + j * m_dim1];
1100 /* L10: */
1101 }
1102 }
1103
1104 /* cleanup odd group of two vectors */
1105
1106 j = *n2 % 4;
1107 if (j >= 2)
1108 {
1109 i__1 = *n1;
1110 for (i__ = 1; i__ <= i__1; ++i__)
1111 {
1112 y[i__] = y[i__] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[i__ + j * m_dim1];
1113 /* L20: */
1114 }
1115 }
1116
1117 /* cleanup odd group of four vectors */
1118
1119 j = *n2 % 8;
1120 if (j >= 4)
1121 {
1122 i__1 = *n1;
1123 for (i__ = 1; i__ <= i__1; ++i__)
1124 {
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];
1126 /* L30: */
1127 }
1128 }
1129
1130 /* cleanup odd group of eight vectors */
1131
1132 j = *n2 % 16;
1133 if (j >= 8)
1134 {
1135 i__1 = *n1;
1136 for (i__ = 1; i__ <= i__1; ++i__)
1137 {
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];
1139 /* L40: */
1140 }
1141 }
1142
1143 /* main loop - groups of sixteen vectors */
1144
1145 jmin = j + 16;
1146 i__1 = *n2;
1147 for (j = jmin; j <= i__1; j += 16)
1148 {
1149 i__2 = *n1;
1150 for (i__ = 1; i__ <= i__2; ++i__)
1151 {
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];
1153 /* L50: */
1154 }
1155 /* L60: */
1156 }
1157 return 0;
1158 } /* dmxpy_ */
1159
ran_(integer * iseed)1160 real_type ran_(integer* iseed)
1161 {
1162 /* System generated locals */
1163 real_type ret_val;
1164
1165 /* Local variables */
1166 static integer it1, it2, it3, it4;
1167
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 */
1173
1174 /* .. Array Arguments .. */
1175 /* .. */
1176
1177 /* Purpose */
1178 /* ======= */
1179
1180 /* DLARAN returns a random double number from a uniform (0,1) */
1181 /* distribution. */
1182
1183 /* Arguments */
1184 /* ========= */
1185
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 */
1189 /* odd. */
1190 /* On exit, the seed is updated. */
1191
1192 /* Further Details */
1193 /* =============== */
1194
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). */
1200
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. */
1204
1205 /* .. Parameters .. */
1206 /* .. */
1207 /* .. Local Scalars .. */
1208 /* .. */
1209 /* .. Intrinsic Functions .. */
1210 /* .. */
1211 /* .. Executable Statements .. */
1212
1213 /* multiply the seed by the multiplier modulo 2**48 */
1214
1215 /* Parameter adjustments */
1216 --iseed;
1217
1218 /* Function Body */
1219 it4 = iseed[4] * 2549;
1220 it3 = it4 / 4096;
1221 it4 -= it3 << 12;
1222 it3 = it3 + iseed[3] * 2549 + iseed[4] * 2508;
1223 it2 = it3 / 4096;
1224 it3 -= it2 << 12;
1225 it2 = it2 + iseed[2] * 2549 + iseed[3] * 2508 + iseed[4] * 322;
1226 it1 = it2 / 4096;
1227 it2 -= it1 << 12;
1228 it1 = it1 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322 + iseed[4] * 494;
1229 it1 %= 4096;
1230
1231 /* return updated seed */
1232
1233 iseed[1] = it1;
1234 iseed[2] = it2;
1235 iseed[3] = it3;
1236 iseed[4] = it4;
1237
1238 /* convert 48-bit integer to a double number in the interval (0,1) */
1239
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;
1241 return ret_val;
1242
1243 /* End of RAN */
1244
1245 } /* ran_ */
1246
1247 /*
1248
1249 Double results:
1250 ~~~~~~~~~~~~~~
1251
1252 norm. resid resid machep x(1) x(n)
1253 6.4915 7.207e-013 2.2204e-016 1 1
1254
1255
1256
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
1261
1262
1263 mpf_class results:
1264 ~~~~~~~~~~~~~~~~~~
1265
1266 norm. resid resid machep x(1) x(n)
1267 3.6575e-05 5.2257e-103 2.8575e-101 1 1
1268
1269
1270
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
1275
1276
1277 number<gmp_float<100> >:
1278 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1279
1280 norm. resid resid machep x(1) x(n)
1281 0.36575e-4 0.52257e-102 0.28575e-100 0.1e1 0.1e1
1282
1283
1284
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
1289
1290 boost::multiprecision::ef::cpp_dec_float_50:
1291 ~~~~~~~~~~~~~~~~~~~~~~~~~
1292
1293 norm. resid resid machep x(1) x(n)
1294 2.551330735e-16 1.275665107e-112 1e-99 1 1
1295
1296
1297
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
1302 */
1303