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 template<>
lexical_cast(const __float128 & f)53 double lexical_cast<double, __float128>(const __float128& f)
54 { return f; }
55
56 }
57
58 #elif defined(TEST_FLOAT128)
59 #include <boost/multiprecision/float128.hpp>
60 typedef boost::multiprecision::float128 real_type;
61 #else
62 typedef double real_type;
63 #endif
64
65 #include <boost/lexical_cast.hpp>
66
67 #ifndef CAST_TO_RT
68 # define CAST_TO_RT(x) x
69 #endif
70
71 extern "C" {
72 #include "f2c.h"
73 integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen),
74 s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
75 e_wsle(void);
76 /* Subroutine */ int s_stop(char *, ftnlen);
77
78 #undef abs
79 #undef dabs
80 #define dabs abs
81 #undef dmin
82 #undef dmax
83 #define dmin min
84 #define dmax max
85
86 }
87 #include <time.h>
88
89 using std::min;
90 using std::max;
91
92 /* Table of constant values */
93
94 static integer c__0 = 0;
95 static real_type c_b7 = CAST_TO_RT(1);
96 static integer c__1 = 1;
97 static integer c__9 = 9;
98
second_(void)99 inline double second_(void)
100 {
101 return ((double)(clock())) / CLOCKS_PER_SEC;
102 }
103
104 int dgefa_(real_type *, integer *, integer *, integer *, integer *), dgesl_(real_type *, integer *, integer *, integer *, real_type *, integer *);
105 int dmxpy_(integer *, real_type *, integer *, integer *, real_type *, real_type *);
106 int matgen_(real_type *, integer *, integer *, real_type *, real_type *);
107 real_type epslon_(real_type *);
108 real_type ran_(integer *);
109 int dscal_(integer *, real_type *, real_type *, integer *);
110 int daxpy_(integer *, real_type *, real_type *, integer *, real_type *, integer *);
111 integer idamax_(integer *, real_type *, integer *);
112 real_type ddot_(integer *, real_type *, integer *, real_type *, integer *);
113 int daxpy_(integer *, real_type *, real_type *, integer *, real_type *, integer *);
114 int dmxpy_(integer *, real_type *, integer *, integer *, real_type *, real_type *);
115
MAIN__()116 extern "C" int MAIN__()
117 {
118 #ifdef TEST_MPF_50
119 std::cout << "Testing number<mpf_float<50> >" << std::endl;
120 #elif defined(TEST_MPFR_50)
121 std::cout << "Testing number<mpf_float<50> >" << std::endl;
122 #elif defined(TEST_GMPXX)
123 std::cout << "Testing mpf_class at 50 decimal degits" << std::endl;
124 mpf_set_default_prec(((50 + 1) * 1000L) / 301L);
125 #elif defined(TEST_MPFRXX)
126 std::cout << "Testing mpfr_class at 50 decimal degits" << std::endl;
127 mpfr_set_default_prec(((50 + 1) * 1000L) / 301L);
128 #elif defined(TEST_CPP_DEC_FLOAT)
129 std::cout << "Testing number<cpp_dec_float<50> >" << std::endl;
130 #elif defined(NATIVE_FLOAT128)
131 std::cout << "Testing __float128" << std::endl;
132 #elif defined(TEST_FLOAT128)
133 std::cout << "Testing number<float128_backend, et_off>" << std::endl;
134 #else
135 std::cout << "Testing double" << std::endl;
136 #endif
137
138
139 /* Format strings */
140 static char fmt_1[] = "(\002 Please send the results of this run to:\002"
141 "//\002 Jack J. Dongarra\002/\002 Computer Science Department\002/"
142 "\002 University of Tennessee\002/\002 Knoxville, Tennessee 37996"
143 "-1300\002//\002 Fax: 615-974-8296\002//\002 Internet: dongarra@c"
144 "s.utk.edu\002/)";
145 static char fmt_40[] = "(\002 norm. resid resid mac"
146 "hep\002,\002 x(1) x(n)\002)";
147 static char fmt_50[] = "(1p5e16.8)";
148 static char fmt_60[] = "(//\002 times are reported for matrices of or"
149 "der \002,i5)";
150 static char fmt_70[] = "(6x,\002factor\002,5x,\002solve\002,6x,\002tota"
151 "l\002,5x,\002mflops\002,7x,\002unit\002,6x,\002ratio\002)";
152 static char fmt_80[] = "(\002 times for array with leading dimension o"
153 "f\002,i4)";
154 static char fmt_110[] = "(6(1pe11.3))";
155
156 /* System generated locals */
157 integer i__1;
158 real_type d__1, d__2, d__3;
159
160 /* Builtin functions */
161
162 /* Local variables */
163 static real_type a[1001000] /* was [1001][1000] */, b[1000];
164 static integer i__, n;
165 static real_type x[1000];
166 static double t1;
167 static integer lda;
168 static double ops;
169 static real_type eps;
170 static integer info;
171 static double time[6], cray, total;
172 static integer ipvt[1000];
173 static real_type resid, norma;
174 static real_type normx;
175 static real_type residn;
176
177 /* Fortran I/O blocks */
178 static cilist io___4 = { 0, 6, 0, fmt_1, 0 };
179 static cilist io___20 = { 0, 6, 0, fmt_40, 0 };
180 static cilist io___21 = { 0, 6, 0, fmt_50, 0 };
181 static cilist io___22 = { 0, 6, 0, fmt_60, 0 };
182 static cilist io___23 = { 0, 6, 0, fmt_70, 0 };
183 static cilist io___24 = { 0, 6, 0, fmt_80, 0 };
184 static cilist io___25 = { 0, 6, 0, fmt_110, 0 };
185 static cilist io___26 = { 0, 6, 0, 0, 0 };
186
187
188 lda = 1001;
189
190 /* this program was updated on 10/12/92 to correct a */
191 /* problem with the random number generator. The previous */
192 /* random number generator had a short period and produced */
193 /* singular matrices occasionally. */
194
195 n = 1000;
196 cray = .056f;
197 s_wsfe(&io___4);
198 e_wsfe();
199 /* Computing 3rd power */
200 d__1 = (real_type) n;
201 /* Computing 2nd power */
202 d__2 = (real_type) n;
203 ops = boost::lexical_cast<double>(real_type(d__1 * (d__1 * d__1) * 2. / 3. + d__2 * d__2 * 2.));
204
205 matgen_(a, &lda, &n, b, &norma);
206
207 /* ****************************************************************** */
208 /* ****************************************************************** */
209 /* you should replace the call to dgefa and dgesl */
210 /* by calls to your linear equation solver. */
211 /* ****************************************************************** */
212 /* ****************************************************************** */
213
214 t1 = second_();
215 dgefa_(a, &lda, &n, ipvt, &info);
216 time[0] = second_() - t1;
217 t1 = second_();
218 dgesl_(a, &lda, &n, ipvt, b, &c__0);
219 time[1] = second_() - t1;
220 total = time[0] + time[1];
221 /* ****************************************************************** */
222 /* ****************************************************************** */
223
224 /* compute a residual to verify results. */
225
226 i__1 = n;
227 for (i__ = 1; i__ <= i__1; ++i__) {
228 x[i__ - 1] = b[i__ - 1];
229 /* L10: */
230 }
231 matgen_(a, &lda, &n, b, &norma);
232 i__1 = n;
233 for (i__ = 1; i__ <= i__1; ++i__) {
234 b[i__ - 1] = -b[i__ - 1];
235 /* L20: */
236 }
237 dmxpy_(&n, b, &n, &lda, x, a);
238 resid = CAST_TO_RT(0);
239 normx = CAST_TO_RT(0);
240 i__1 = n;
241 for (i__ = 1; i__ <= i__1; ++i__) {
242 /* Computing MAX */
243 d__2 = resid, d__3 = (d__1 = b[i__ - 1], abs(d__1));
244 resid = (max)(d__2,d__3);
245 /* Computing MAX */
246 d__2 = normx, d__3 = (d__1 = x[i__ - 1], abs(d__1));
247 normx = (max)(d__2,d__3);
248 /* L30: */
249 }
250 eps = epslon_(&c_b7);
251 residn = resid / (n * norma * normx * eps);
252 s_wsfe(&io___20);
253 e_wsfe();
254 s_wsfe(&io___21);
255 /*
256 do_fio(&c__1, (char *)&residn, (ftnlen)sizeof(real_type));
257 do_fio(&c__1, (char *)&resid, (ftnlen)sizeof(real_type));
258 do_fio(&c__1, (char *)&eps, (ftnlen)sizeof(real_type));
259 do_fio(&c__1, (char *)&x[0], (ftnlen)sizeof(real_type));
260 do_fio(&c__1, (char *)&x[n - 1], (ftnlen)sizeof(real_type));
261 */
262 std::cout << std::setw(12) << std::setprecision(5) << residn << " " << resid << " " << eps << " " << x[0] << " " << x[n-1] << std::endl;
263 e_wsfe();
264
265 s_wsfe(&io___22);
266 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
267 e_wsfe();
268 s_wsfe(&io___23);
269 e_wsfe();
270
271 time[2] = total;
272 time[3] = ops / (total * 1e6);
273 time[4] = 2. / time[3];
274 time[5] = total / cray;
275 s_wsfe(&io___24);
276 do_fio(&c__1, (char *)&lda, (ftnlen)sizeof(integer));
277 e_wsfe();
278 s_wsfe(&io___25);
279 for (i__ = 1; i__ <= 6; ++i__) {
280 // do_fio(&c__1, (char *)&time[i__ - 1], (ftnlen)sizeof(real_type));
281 std::cout << std::setw(12) << std::setprecision(5) << time[i__ - 1];
282 }
283 e_wsfe();
284 s_wsle(&io___26);
285 do_lio(&c__9, &c__1, " end of tests -- this version dated 10/12/92", (
286 ftnlen)44);
287 e_wsle();
288
289 s_stop("", (ftnlen)0);
290 return 0;
291 } /* MAIN__ */
292
matgen_(real_type * a,integer * lda,integer * n,real_type * b,real_type * norma)293 /* Subroutine */ int matgen_(real_type *a, integer *lda, integer *n,
294 real_type *b, real_type *norma)
295 {
296 /* System generated locals */
297 integer a_dim1, a_offset, i__1, i__2;
298 real_type d__1, d__2;
299
300 /* Local variables */
301 static integer i__, j;
302 static integer init[4];
303
304
305 /* Parameter adjustments */
306 a_dim1 = *lda;
307 a_offset = 1 + a_dim1;
308 a -= a_offset;
309 --b;
310
311 /* Function Body */
312 init[0] = 1;
313 init[1] = 2;
314 init[2] = 3;
315 init[3] = 1325;
316 *norma = CAST_TO_RT(0);
317 i__1 = *n;
318 for (j = 1; j <= i__1; ++j) {
319 i__2 = *n;
320 for (i__ = 1; i__ <= i__2; ++i__) {
321 a[i__ + j * a_dim1] = ran_(init) - .5f;
322 /* Computing MAX */
323 d__2 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
324 *norma = (max)(d__2,*norma);
325 /* L20: */
326 }
327 /* L30: */
328 }
329 i__1 = *n;
330 for (i__ = 1; i__ <= i__1; ++i__) {
331 b[i__] = CAST_TO_RT(0);
332 /* L35: */
333 }
334 i__1 = *n;
335 for (j = 1; j <= i__1; ++j) {
336 i__2 = *n;
337 for (i__ = 1; i__ <= i__2; ++i__) {
338 b[i__] += a[i__ + j * a_dim1];
339 /* L40: */
340 }
341 /* L50: */
342 }
343 return 0;
344 } /* matgen_ */
345
dgefa_(real_type * a,integer * lda,integer * n,integer * ipvt,integer * info)346 /* Subroutine */ int dgefa_(real_type *a, integer *lda, integer *n, integer *
347 ipvt, integer *info)
348 {
349 /* System generated locals */
350 integer a_dim1, a_offset, i__1, i__2, i__3;
351
352 /* Local variables */
353 static integer j, k, l;
354 static real_type t;
355 static integer kp1, nm1;
356
357
358 /* dgefa factors a double precision matrix by gaussian elimination. */
359
360 /* dgefa is usually called by dgeco, but it can be called */
361 /* directly with a saving in time if rcond is not needed. */
362 /* (time for dgeco) = (1 + 9/n)*(time for dgefa) . */
363
364 /* on entry */
365
366 /* a double precision(lda, n) */
367 /* the matrix to be factored. */
368
369 /* lda integer */
370 /* the leading dimension of the array a . */
371
372 /* n integer */
373 /* the order of the matrix a . */
374
375 /* on return */
376
377 /* a an upper triangular matrix and the multipliers */
378 /* which were used to obtain it. */
379 /* the factorization can be written a = l*u where */
380 /* l is a product of permutation and unit lower */
381 /* triangular matrices and u is upper triangular. */
382
383 /* ipvt integer(n) */
384 /* an integer vector of pivot indices. */
385
386 /* info integer */
387 /* = 0 normal value. */
388 /* = k if u(k,k) .eq. 0.0 . this is not an error */
389 /* condition for this subroutine, but it does */
390 /* indicate that dgesl or dgedi will divide by zero */
391 /* if called. use rcond in dgeco for a reliable */
392 /* indication of singularity. */
393
394 /* linpack. this version dated 08/14/78 . */
395 /* cleve moler, university of new mexico, argonne national lab. */
396
397 /* subroutines and functions */
398
399 /* blas daxpy,dscal,idamax */
400
401 /* internal variables */
402
403
404
405 /* gaussian elimination with partial pivoting */
406
407 /* Parameter adjustments */
408 a_dim1 = *lda;
409 a_offset = 1 + a_dim1;
410 a -= a_offset;
411 --ipvt;
412
413 /* Function Body */
414 *info = 0;
415 nm1 = *n - 1;
416 if (nm1 < 1) {
417 goto L70;
418 }
419 i__1 = nm1;
420 for (k = 1; k <= i__1; ++k) {
421 kp1 = k + 1;
422
423 /* find l = pivot index */
424
425 i__2 = *n - k + 1;
426 l = idamax_(&i__2, &a[k + k * a_dim1], &c__1) + k - 1;
427 ipvt[k] = l;
428
429 /* zero pivot implies this column already triangularized */
430
431 if (a[l + k * a_dim1] == 0.) {
432 goto L40;
433 }
434
435 /* interchange if necessary */
436
437 if (l == k) {
438 goto L10;
439 }
440 t = a[l + k * a_dim1];
441 a[l + k * a_dim1] = a[k + k * a_dim1];
442 a[k + k * a_dim1] = t;
443 L10:
444
445 /* compute multipliers */
446
447 t = -1. / a[k + k * a_dim1];
448 i__2 = *n - k;
449 dscal_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1);
450
451 /* row elimination with column indexing */
452
453 i__2 = *n;
454 for (j = kp1; j <= i__2; ++j) {
455 t = a[l + j * a_dim1];
456 if (l == k) {
457 goto L20;
458 }
459 a[l + j * a_dim1] = a[k + j * a_dim1];
460 a[k + j * a_dim1] = t;
461 L20:
462 i__3 = *n - k;
463 daxpy_(&i__3, &t, &a[k + 1 + k * a_dim1], &c__1, &a[k + 1 + j *
464 a_dim1], &c__1);
465 /* L30: */
466 }
467 goto L50;
468 L40:
469 *info = k;
470 L50:
471 /* L60: */
472 ;
473 }
474 L70:
475 ipvt[*n] = *n;
476 if (a[*n + *n * a_dim1] == 0.) {
477 *info = *n;
478 }
479 return 0;
480 } /* dgefa_ */
481
dgesl_(real_type * a,integer * lda,integer * n,integer * ipvt,real_type * b,integer * job)482 /* Subroutine */ int dgesl_(real_type *a, integer *lda, integer *n, integer *
483 ipvt, real_type *b, integer *job)
484 {
485 /* System generated locals */
486 integer a_dim1, a_offset, i__1, i__2;
487
488 /* Local variables */
489 static integer k, l;
490 static real_type t;
491 static integer kb, nm1;
492
493
494 /* dgesl solves the double precision system */
495 /* a * x = b or trans(a) * x = b */
496 /* using the factors computed by dgeco or dgefa. */
497
498 /* on entry */
499
500 /* a double precision(lda, n) */
501 /* the output from dgeco or dgefa. */
502
503 /* lda integer */
504 /* the leading dimension of the array a . */
505
506 /* n integer */
507 /* the order of the matrix a . */
508
509 /* ipvt integer(n) */
510 /* the pivot vector from dgeco or dgefa. */
511
512 /* b double precision(n) */
513 /* the right hand side vector. */
514
515 /* job integer */
516 /* = 0 to solve a*x = b , */
517 /* = nonzero to solve trans(a)*x = b where */
518 /* trans(a) is the transpose. */
519
520 /* on return */
521
522 /* b the solution vector x . */
523
524 /* error condition */
525
526 /* a division by zero will occur if the input factor contains a */
527 /* zero on the diagonal. technically this indicates singularity */
528 /* but it is often caused by improper arguments or improper */
529 /* setting of lda . it will not occur if the subroutines are */
530 /* called correctly and if dgeco has set rcond .gt. 0.0 */
531 /* or dgefa has set info .eq. 0 . */
532
533 /* to compute inverse(a) * c where c is a matrix */
534 /* with p columns */
535 /* call dgeco(a,lda,n,ipvt,rcond,z) */
536 /* if (rcond is too small) go to ... */
537 /* do 10 j = 1, p */
538 /* call dgesl(a,lda,n,ipvt,c(1,j),0) */
539 /* 10 continue */
540
541 /* linpack. this version dated 08/14/78 . */
542 /* cleve moler, university of new mexico, argonne national lab. */
543
544 /* subroutines and functions */
545
546 /* blas daxpy,ddot */
547
548 /* internal variables */
549
550
551 /* Parameter adjustments */
552 a_dim1 = *lda;
553 a_offset = 1 + a_dim1;
554 a -= a_offset;
555 --ipvt;
556 --b;
557
558 /* Function Body */
559 nm1 = *n - 1;
560 if (*job != 0) {
561 goto L50;
562 }
563
564 /* job = 0 , solve a * x = b */
565 /* first solve l*y = b */
566
567 if (nm1 < 1) {
568 goto L30;
569 }
570 i__1 = nm1;
571 for (k = 1; k <= i__1; ++k) {
572 l = ipvt[k];
573 t = b[l];
574 if (l == k) {
575 goto L10;
576 }
577 b[l] = b[k];
578 b[k] = t;
579 L10:
580 i__2 = *n - k;
581 daxpy_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
582 /* L20: */
583 }
584 L30:
585
586 /* now solve u*x = y */
587
588 i__1 = *n;
589 for (kb = 1; kb <= i__1; ++kb) {
590 k = *n + 1 - kb;
591 b[k] /= a[k + k * a_dim1];
592 t = -b[k];
593 i__2 = k - 1;
594 daxpy_(&i__2, &t, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
595 /* L40: */
596 }
597 goto L100;
598 L50:
599
600 /* job = nonzero, solve trans(a) * x = b */
601 /* first solve trans(u)*y = b */
602
603 i__1 = *n;
604 for (k = 1; k <= i__1; ++k) {
605 i__2 = k - 1;
606 t = ddot_(&i__2, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
607 b[k] = (b[k] - t) / a[k + k * a_dim1];
608 /* L60: */
609 }
610
611 /* now solve trans(l)*x = y */
612
613 if (nm1 < 1) {
614 goto L90;
615 }
616 i__1 = nm1;
617 for (kb = 1; kb <= i__1; ++kb) {
618 k = *n - kb;
619 i__2 = *n - k;
620 b[k] += ddot_(&i__2, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
621 l = ipvt[k];
622 if (l == k) {
623 goto L70;
624 }
625 t = b[l];
626 b[l] = b[k];
627 b[k] = t;
628 L70:
629 /* L80: */
630 ;
631 }
632 L90:
633 L100:
634 return 0;
635 } /* dgesl_ */
636
daxpy_(integer * n,real_type * da,real_type * dx,integer * incx,real_type * dy,integer * incy)637 /* Subroutine */ int daxpy_(integer *n, real_type *da, real_type *dx,
638 integer *incx, real_type *dy, integer *incy)
639 {
640 /* System generated locals */
641 integer i__1;
642
643 /* Local variables */
644 static integer i__, m, ix, iy, mp1;
645
646
647 /* constant times a vector plus a vector. */
648 /* uses unrolled loops for increments equal to one. */
649 /* jack dongarra, linpack, 3/11/78. */
650
651
652 /* Parameter adjustments */
653 --dy;
654 --dx;
655
656 /* Function Body */
657 if (*n <= 0) {
658 return 0;
659 }
660 if (*da == 0.) {
661 return 0;
662 }
663 if (*incx == 1 && *incy == 1) {
664 goto L20;
665 }
666
667 /* code for unequal increments or equal increments */
668 /* not equal to 1 */
669
670 ix = 1;
671 iy = 1;
672 if (*incx < 0) {
673 ix = (-(*n) + 1) * *incx + 1;
674 }
675 if (*incy < 0) {
676 iy = (-(*n) + 1) * *incy + 1;
677 }
678 i__1 = *n;
679 for (i__ = 1; i__ <= i__1; ++i__) {
680 dy[iy] += *da * dx[ix];
681 ix += *incx;
682 iy += *incy;
683 /* L10: */
684 }
685 return 0;
686
687 /* code for both increments equal to 1 */
688
689
690 /* clean-up loop */
691
692 L20:
693 m = *n % 4;
694 if (m == 0) {
695 goto L40;
696 }
697 i__1 = m;
698 for (i__ = 1; i__ <= i__1; ++i__) {
699 dy[i__] += *da * dx[i__];
700 /* L30: */
701 }
702 if (*n < 4) {
703 return 0;
704 }
705 L40:
706 mp1 = m + 1;
707 i__1 = *n;
708 for (i__ = mp1; i__ <= i__1; i__ += 4) {
709 dy[i__] += *da * dx[i__];
710 dy[i__ + 1] += *da * dx[i__ + 1];
711 dy[i__ + 2] += *da * dx[i__ + 2];
712 dy[i__ + 3] += *da * dx[i__ + 3];
713 /* L50: */
714 }
715 return 0;
716 } /* daxpy_ */
717
ddot_(integer * n,real_type * dx,integer * incx,real_type * dy,integer * incy)718 real_type ddot_(integer *n, real_type *dx, integer *incx, real_type *dy,
719 integer *incy)
720 {
721 /* System generated locals */
722 integer i__1;
723 real_type ret_val;
724
725 /* Local variables */
726 static integer i__, m, ix, iy, mp1;
727 static real_type dtemp;
728
729
730 /* forms the dot product of two vectors. */
731 /* uses unrolled loops for increments equal to one. */
732 /* jack dongarra, linpack, 3/11/78. */
733
734
735 /* Parameter adjustments */
736 --dy;
737 --dx;
738
739 /* Function Body */
740 ret_val = CAST_TO_RT(0);
741 dtemp = CAST_TO_RT(0);
742 if (*n <= 0) {
743 return ret_val;
744 }
745 if (*incx == 1 && *incy == 1) {
746 goto L20;
747 }
748
749 /* code for unequal increments or equal increments */
750 /* not equal to 1 */
751
752 ix = 1;
753 iy = 1;
754 if (*incx < 0) {
755 ix = (-(*n) + 1) * *incx + 1;
756 }
757 if (*incy < 0) {
758 iy = (-(*n) + 1) * *incy + 1;
759 }
760 i__1 = *n;
761 for (i__ = 1; i__ <= i__1; ++i__) {
762 dtemp += dx[ix] * dy[iy];
763 ix += *incx;
764 iy += *incy;
765 /* L10: */
766 }
767 ret_val = dtemp;
768 return ret_val;
769
770 /* code for both increments equal to 1 */
771
772
773 /* clean-up loop */
774
775 L20:
776 m = *n % 5;
777 if (m == 0) {
778 goto L40;
779 }
780 i__1 = m;
781 for (i__ = 1; i__ <= i__1; ++i__) {
782 dtemp += dx[i__] * dy[i__];
783 /* L30: */
784 }
785 if (*n < 5) {
786 goto L60;
787 }
788 L40:
789 mp1 = m + 1;
790 i__1 = *n;
791 for (i__ = mp1; i__ <= i__1; i__ += 5) {
792 dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
793 i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
794 4] * dy[i__ + 4];
795 /* L50: */
796 }
797 L60:
798 ret_val = dtemp;
799 return ret_val;
800 } /* ddot_ */
801
dscal_(integer * n,real_type * da,real_type * dx,integer * incx)802 /* Subroutine */ int dscal_(integer *n, real_type *da, real_type *dx,
803 integer *incx)
804 {
805 /* System generated locals */
806 integer i__1, i__2;
807
808 /* Local variables */
809 static integer i__, m, mp1, nincx;
810
811
812 /* scales a vector by a constant. */
813 /* uses unrolled loops for increment equal to one. */
814 /* jack dongarra, linpack, 3/11/78. */
815
816
817 /* Parameter adjustments */
818 --dx;
819
820 /* Function Body */
821 if (*n <= 0) {
822 return 0;
823 }
824 if (*incx == 1) {
825 goto L20;
826 }
827
828 /* code for increment not equal to 1 */
829
830 nincx = *n * *incx;
831 i__1 = nincx;
832 i__2 = *incx;
833 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
834 dx[i__] = *da * dx[i__];
835 /* L10: */
836 }
837 return 0;
838
839 /* code for increment equal to 1 */
840
841
842 /* clean-up loop */
843
844 L20:
845 m = *n % 5;
846 if (m == 0) {
847 goto L40;
848 }
849 i__2 = m;
850 for (i__ = 1; i__ <= i__2; ++i__) {
851 dx[i__] = *da * dx[i__];
852 /* L30: */
853 }
854 if (*n < 5) {
855 return 0;
856 }
857 L40:
858 mp1 = m + 1;
859 i__2 = *n;
860 for (i__ = mp1; i__ <= i__2; i__ += 5) {
861 dx[i__] = *da * dx[i__];
862 dx[i__ + 1] = *da * dx[i__ + 1];
863 dx[i__ + 2] = *da * dx[i__ + 2];
864 dx[i__ + 3] = *da * dx[i__ + 3];
865 dx[i__ + 4] = *da * dx[i__ + 4];
866 /* L50: */
867 }
868 return 0;
869 } /* dscal_ */
870
idamax_(integer * n,real_type * dx,integer * incx)871 integer idamax_(integer *n, real_type *dx, integer *incx)
872 {
873 /* System generated locals */
874 integer ret_val, i__1;
875 real_type d__1;
876
877 /* Local variables */
878 static integer i__, ix;
879 static real_type dmax__;
880
881
882 /* finds the index of element having max. dabsolute value. */
883 /* jack dongarra, linpack, 3/11/78. */
884
885
886 /* Parameter adjustments */
887 --dx;
888
889 /* Function Body */
890 ret_val = 0;
891 if (*n < 1) {
892 return ret_val;
893 }
894 ret_val = 1;
895 if (*n == 1) {
896 return ret_val;
897 }
898 if (*incx == 1) {
899 goto L20;
900 }
901
902 /* code for increment not equal to 1 */
903
904 ix = 1;
905 dmax__ = abs(dx[1]);
906 ix += *incx;
907 i__1 = *n;
908 for (i__ = 2; i__ <= i__1; ++i__) {
909 if ((d__1 = dx[ix], abs(d__1)) <= dmax__) {
910 goto L5;
911 }
912 ret_val = i__;
913 dmax__ = (d__1 = dx[ix], abs(d__1));
914 L5:
915 ix += *incx;
916 /* L10: */
917 }
918 return ret_val;
919
920 /* code for increment equal to 1 */
921
922 L20:
923 dmax__ = abs(dx[1]);
924 i__1 = *n;
925 for (i__ = 2; i__ <= i__1; ++i__) {
926 if ((d__1 = dx[i__], abs(d__1)) <= dmax__) {
927 goto L30;
928 }
929 ret_val = i__;
930 dmax__ = (d__1 = dx[i__], abs(d__1));
931 L30:
932 ;
933 }
934 return ret_val;
935 } /* idamax_ */
936
epslon_(real_type * x)937 real_type epslon_(real_type *x)
938 {
939 #if defined(TEST_MPF_100) || defined(TEST_MPFR_100) || defined(TEST_GMPXX) || defined(TEST_MPFRXX)
940 return std::ldexp(1.0, 1 - ((100 + 1) * 1000L) / 301L);
941 #elif defined(TEST_CPP_DEC_FLOAT_BN)
942 return std::pow(10.0, 1-std::numeric_limits<efx::cpp_dec_float_50>::digits10);
943 #elif defined(NATIVE_FLOAT128)
944 return FLT128_EPSILON;
945 #else
946 return CAST_TO_RT(std::numeric_limits<real_type>::epsilon());
947 #endif
948 } /* epslon_ */
949
mm_(real_type * a,integer * lda,integer * n1,integer * n3,real_type * b,integer * ldb,integer * n2,real_type * c__,integer * ldc)950 /* Subroutine */ int mm_(real_type *a, integer *lda, integer *n1, integer *
951 n3, real_type *b, integer *ldb, integer *n2, real_type *c__,
952 integer *ldc)
953 {
954 /* System generated locals */
955 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2;
956
957 /* Local variables */
958 static integer i__, j;
959
960
961 /* purpose: */
962 /* multiply matrix b times matrix c and store the result in matrix a. */
963
964 /* parameters: */
965
966 /* a double precision(lda,n3), matrix of n1 rows and n3 columns */
967
968 /* lda integer, leading dimension of array a */
969
970 /* n1 integer, number of rows in matrices a and b */
971
972 /* n3 integer, number of columns in matrices a and c */
973
974 /* b double precision(ldb,n2), matrix of n1 rows and n2 columns */
975
976 /* ldb integer, leading dimension of array b */
977
978 /* n2 integer, number of columns in matrix b, and number of rows in */
979 /* matrix c */
980
981 /* c double precision(ldc,n3), matrix of n2 rows and n3 columns */
982
983 /* ldc integer, leading dimension of array c */
984
985 /* ---------------------------------------------------------------------- */
986
987 /* Parameter adjustments */
988 a_dim1 = *lda;
989 a_offset = 1 + a_dim1;
990 a -= a_offset;
991 b_dim1 = *ldb;
992 b_offset = 1 + b_dim1;
993 b -= b_offset;
994 c_dim1 = *ldc;
995 c_offset = 1 + c_dim1;
996 c__ -= c_offset;
997
998 /* Function Body */
999 i__1 = *n3;
1000 for (j = 1; j <= i__1; ++j) {
1001 i__2 = *n1;
1002 for (i__ = 1; i__ <= i__2; ++i__) {
1003 a[i__ + j * a_dim1] = CAST_TO_RT(0);
1004 /* L10: */
1005 }
1006 dmxpy_(n2, &a[j * a_dim1 + 1], n1, ldb, &c__[j * c_dim1 + 1], &b[
1007 b_offset]);
1008 /* L20: */
1009 }
1010
1011 return 0;
1012 } /* mm_ */
1013
dmxpy_(integer * n1,real_type * y,integer * n2,integer * ldm,real_type * x,real_type * m)1014 /* Subroutine */ int dmxpy_(integer *n1, real_type *y, integer *n2, integer *
1015 ldm, real_type *x, real_type *m)
1016 {
1017 /* System generated locals */
1018 integer m_dim1, m_offset, i__1, i__2;
1019
1020 /* Local variables */
1021 static integer i__, j, jmin;
1022
1023
1024 /* purpose: */
1025 /* multiply matrix m times vector x and add the result to vector y. */
1026
1027 /* parameters: */
1028
1029 /* n1 integer, number of elements in vector y, and number of rows in */
1030 /* matrix m */
1031
1032 /* y double precision(n1), vector of length n1 to which is added */
1033 /* the product m*x */
1034
1035 /* n2 integer, number of elements in vector x, and number of columns */
1036 /* in matrix m */
1037
1038 /* ldm integer, leading dimension of array m */
1039
1040 /* x double precision(n2), vector of length n2 */
1041
1042 /* m double precision(ldm,n2), matrix of n1 rows and n2 columns */
1043
1044 /* ---------------------------------------------------------------------- */
1045
1046 /* cleanup odd vector */
1047
1048 /* Parameter adjustments */
1049 --y;
1050 m_dim1 = *ldm;
1051 m_offset = 1 + m_dim1;
1052 m -= m_offset;
1053 --x;
1054
1055 /* Function Body */
1056 j = *n2 % 2;
1057 if (j >= 1) {
1058 i__1 = *n1;
1059 for (i__ = 1; i__ <= i__1; ++i__) {
1060 y[i__] += x[j] * m[i__ + j * m_dim1];
1061 /* L10: */
1062 }
1063 }
1064
1065 /* cleanup odd group of two vectors */
1066
1067 j = *n2 % 4;
1068 if (j >= 2) {
1069 i__1 = *n1;
1070 for (i__ = 1; i__ <= i__1; ++i__) {
1071 y[i__] = y[i__] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[
1072 i__ + j * m_dim1];
1073 /* L20: */
1074 }
1075 }
1076
1077 /* cleanup odd group of four vectors */
1078
1079 j = *n2 % 8;
1080 if (j >= 4) {
1081 i__1 = *n1;
1082 for (i__ = 1; i__ <= i__1; ++i__) {
1083 y[i__] = y[i__] + x[j - 3] * m[i__ + (j - 3) * m_dim1] + x[j - 2]
1084 * m[i__ + (j - 2) * m_dim1] + x[j - 1] * m[i__ + (j - 1) *
1085 m_dim1] + x[j] * m[i__ + j * m_dim1];
1086 /* L30: */
1087 }
1088 }
1089
1090 /* cleanup odd group of eight vectors */
1091
1092 j = *n2 % 16;
1093 if (j >= 8) {
1094 i__1 = *n1;
1095 for (i__ = 1; i__ <= i__1; ++i__) {
1096 y[i__] = y[i__] + x[j - 7] * m[i__ + (j - 7) * m_dim1] + x[j - 6]
1097 * m[i__ + (j - 6) * m_dim1] + x[j - 5] * m[i__ + (j - 5) *
1098 m_dim1] + x[j - 4] * m[i__ + (j - 4) * m_dim1] + x[j - 3]
1099 * m[i__ + (j - 3) * m_dim1] + x[j - 2] * m[i__ + (j - 2)
1100 * m_dim1] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] *
1101 m[i__ + j * m_dim1];
1102 /* L40: */
1103 }
1104 }
1105
1106 /* main loop - groups of sixteen vectors */
1107
1108 jmin = j + 16;
1109 i__1 = *n2;
1110 for (j = jmin; j <= i__1; j += 16) {
1111 i__2 = *n1;
1112 for (i__ = 1; i__ <= i__2; ++i__) {
1113 y[i__] = y[i__] + x[j - 15] * m[i__ + (j - 15) * m_dim1] + x[j -
1114 14] * m[i__ + (j - 14) * m_dim1] + x[j - 13] * m[i__ + (j
1115 - 13) * m_dim1] + x[j - 12] * m[i__ + (j - 12) * m_dim1]
1116 + x[j - 11] * m[i__ + (j - 11) * m_dim1] + x[j - 10] * m[
1117 i__ + (j - 10) * m_dim1] + x[j - 9] * m[i__ + (j - 9) *
1118 m_dim1] + x[j - 8] * m[i__ + (j - 8) * m_dim1] + x[j - 7]
1119 * m[i__ + (j - 7) * m_dim1] + x[j - 6] * m[i__ + (j - 6) *
1120 m_dim1] + x[j - 5] * m[i__ + (j - 5) * m_dim1] + x[j - 4]
1121 * m[i__ + (j - 4) * m_dim1] + x[j - 3] * m[i__ + (j - 3)
1122 * m_dim1] + x[j - 2] * m[i__ + (j - 2) * m_dim1] + x[j -
1123 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[i__ + j *
1124 m_dim1];
1125 /* L50: */
1126 }
1127 /* L60: */
1128 }
1129 return 0;
1130 } /* dmxpy_ */
1131
ran_(integer * iseed)1132 real_type ran_(integer *iseed)
1133 {
1134 /* System generated locals */
1135 real_type ret_val;
1136
1137 /* Local variables */
1138 static integer it1, it2, it3, it4;
1139
1140
1141 /* modified from the LAPACK auxiliary routine 10/12/92 JD */
1142 /* -- LAPACK auxiliary routine (version 1.0) -- */
1143 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
1144 /* Courant Institute, Argonne National Lab, and Rice University */
1145 /* February 29, 1992 */
1146
1147 /* .. Array Arguments .. */
1148 /* .. */
1149
1150 /* Purpose */
1151 /* ======= */
1152
1153 /* DLARAN returns a random double number from a uniform (0,1) */
1154 /* distribution. */
1155
1156 /* Arguments */
1157 /* ========= */
1158
1159 /* ISEED (input/output) INTEGER array, dimension (4) */
1160 /* On entry, the seed of the random number generator; the array */
1161 /* elements must be between 0 and 4095, and ISEED(4) must be */
1162 /* odd. */
1163 /* On exit, the seed is updated. */
1164
1165 /* Further Details */
1166 /* =============== */
1167
1168 /* This routine uses a multiplicative congruential method with modulus */
1169 /* 2**48 and multiplier 33952834046453 (see G.S.Fishman, */
1170 /* 'Multiplicative congruential random number generators with modulus */
1171 /* 2**b: an exhaustive analysis for b = 32 and a partial analysis for */
1172 /* b = 48', Math. Comp. 189, pp 331-344, 1990). */
1173
1174 /* 48-bit integers are stored in 4 integer array elements with 12 bits */
1175 /* per element. Hence the routine is portable across machines with */
1176 /* integers of 32 bits or more. */
1177
1178 /* .. Parameters .. */
1179 /* .. */
1180 /* .. Local Scalars .. */
1181 /* .. */
1182 /* .. Intrinsic Functions .. */
1183 /* .. */
1184 /* .. Executable Statements .. */
1185
1186 /* multiply the seed by the multiplier modulo 2**48 */
1187
1188 /* Parameter adjustments */
1189 --iseed;
1190
1191 /* Function Body */
1192 it4 = iseed[4] * 2549;
1193 it3 = it4 / 4096;
1194 it4 -= it3 << 12;
1195 it3 = it3 + iseed[3] * 2549 + iseed[4] * 2508;
1196 it2 = it3 / 4096;
1197 it3 -= it2 << 12;
1198 it2 = it2 + iseed[2] * 2549 + iseed[3] * 2508 + iseed[4] * 322;
1199 it1 = it2 / 4096;
1200 it2 -= it1 << 12;
1201 it1 = it1 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322 + iseed[4]
1202 * 494;
1203 it1 %= 4096;
1204
1205 /* return updated seed */
1206
1207 iseed[1] = it1;
1208 iseed[2] = it2;
1209 iseed[3] = it3;
1210 iseed[4] = it4;
1211
1212 /* convert 48-bit integer to a double number in the interval (0,1) */
1213
1214 ret_val = ((real_type) it1 + ((real_type) it2 + ((real_type) it3 + (
1215 real_type) it4 * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4)
1216 * 2.44140625e-4;
1217 return ret_val;
1218
1219 /* End of RAN */
1220
1221 } /* ran_ */
1222
1223 /*
1224
1225 Double results:
1226 ~~~~~~~~~~~~~~
1227
1228 norm. resid resid machep x(1) x(n)
1229 6.4915 7.207e-013 2.2204e-016 1 1
1230
1231
1232
1233 times are reported for matrices of order 1000
1234 factor solve total mflops unit ratio
1235 times for array with leading dimension of1001
1236 1.443 0.003 1.446 462.43 0.004325 25.821
1237
1238
1239 mpf_class results:
1240 ~~~~~~~~~~~~~~~~~~
1241
1242 norm. resid resid machep x(1) x(n)
1243 3.6575e-05 5.2257e-103 2.8575e-101 1 1
1244
1245
1246
1247 times are reported for matrices of order 1000
1248 factor solve total mflops unit ratio
1249 times for array with leading dimension of1001
1250 266.45 0.798 267.24 2.5021 0.79933 4772.2
1251
1252
1253 number<gmp_float<100> >:
1254 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1255
1256 norm. resid resid machep x(1) x(n)
1257 0.36575e-4 0.52257e-102 0.28575e-100 0.1e1 0.1e1
1258
1259
1260
1261 times are reported for matrices of order 1000
1262 factor solve total mflops unit ratio
1263 times for array with leading dimension of1001
1264 279.96 0.84 280.8 2.3813 0.83988 5014.3
1265
1266 boost::multiprecision::ef::cpp_dec_float_50:
1267 ~~~~~~~~~~~~~~~~~~~~~~~~~
1268
1269 norm. resid resid machep x(1) x(n)
1270 2.551330735e-16 1.275665107e-112 1e-99 1 1
1271
1272
1273
1274 times are reported for matrices of order 1000
1275 factor solve total mflops unit ratio
1276 times for array with leading dimension of1001
1277 363.89 1.074 364.97 1.8321 1.0916 6517.3
1278 */
1279