1 # include <stdlib.h>
2 # include <stdio.h>
3 # include <math.h>
4 # include <time.h>
5
6 double cpu_time ( );
7 void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy );
8 double ddot ( int n, double dx[], int incx, double dy[], int incy );
9 int dgefa ( double a[], int lda, int n, int ipvt[] );
10 void dgesl ( double a[], int lda, int n, int ipvt[], double b[], int job );
11 void dscal ( int n, double sa, double x[], int incx );
12 int idamax ( int n, double dx[], int incx );
13 double r8_abs ( double x );
14 double r8_epsilon ( );
15 double r8_max ( double x, double y );
16 double r8_random ( int iseed[4] );
17 double *r8mat_gen ( int lda, int n );
18 void timestamp ( );
19
20 /******************************************************************************/
21
main(int argc,char ** argv)22 int main(int argc, char **argv) {
23
24 /******************************************************************************/
25 /*
26 Purpose:
27
28 MAIN is the main program for LINPACK_BENCH.
29
30 Discussion:
31
32 LINPACK_BENCH drives the double precision LINPACK benchmark program.
33
34 Modified:
35
36 25 July 2008
37
38 Parameters:
39
40 N is the problem size.
41 */
42 # define N 2000
43 # define LDA ( N + 1 )
44
45 double *a;
46 double a_max;
47 double *b;
48 double b_max;
49 double cray = 0.056;
50 double eps;
51 int i;
52 int info;
53 int *ipvt;
54 int j;
55 int job;
56 double ops;
57 double *resid;
58 double resid_max;
59 double residn;
60 double *rhs;
61 double t1;
62 double t2;
63 double time[6];
64 double total;
65 double *x;
66
67 int arg = argc > 1 ? argv[1][0] - '0' : 3;
68 if (arg == 0) return 0;
69
70 timestamp ( );
71 printf ( "\n" );
72 printf ( "LINPACK_BENCH\n" );
73 printf ( " C version\n" );
74 printf ( "\n" );
75 printf ( " The LINPACK benchmark.\n" );
76 printf ( " Language: C\n" );
77 printf ( " Datatype: Double precision real\n" );
78 printf ( " Matrix order N = %d\n", N );
79 printf ( " Leading matrix dimension LDA = %d\n", LDA );
80
81 ops = ( double ) ( 2.0 * N * N * N ) / 3.0 + 2.0 * ( double ) ( N * N );
82 /*
83 Allocate space for arrays.
84 */
85 a = r8mat_gen ( LDA, N );
86 b = ( double * ) malloc ( N * sizeof ( double ) );
87 ipvt = ( int * ) malloc ( N * sizeof ( int ) );
88 resid = ( double * ) malloc ( N * sizeof ( double ) );
89 rhs = ( double * ) malloc ( N * sizeof ( double ) );
90 x = ( double * ) malloc ( N * sizeof ( double ) );
91
92 a_max = 0.0;
93 for ( j = 0; j < N; j++ )
94 {
95 for ( i = 0; i < N; i++ )
96 {
97 a_max = r8_max ( a_max, a[i+j*LDA] );
98 }
99 }
100
101 for ( i = 0; i < N; i++ )
102 {
103 x[i] = 1.0;
104 }
105
106 for ( i = 0; i < N; i++ )
107 {
108 b[i] = 0.0;
109 for ( j = 0; j < N; j++ )
110 {
111 b[i] = b[i] + a[i+j*LDA] * x[j];
112 }
113 }
114 t1 = cpu_time ( );
115
116 info = dgefa ( a, LDA, N, ipvt );
117
118 if ( info != 0 )
119 {
120 printf ( "\n" );
121 printf ( "LINPACK_BENCH - Fatal error!\n" );
122 printf ( " The matrix A is apparently singular.\n" );
123 printf ( " Abnormal end of execution.\n" );
124 return 1;
125 }
126
127 t2 = cpu_time ( );
128 time[0] = t2 - t1;
129
130 t1 = cpu_time ( );
131
132 job = 0;
133 dgesl ( a, LDA, N, ipvt, b, job );
134
135 t2 = cpu_time ( );
136 time[1] = t2 - t1;
137
138 total = time[0] + time[1];
139
140 free ( a );
141 /*
142 Compute a residual to verify results.
143 */
144 a = r8mat_gen ( LDA, N );
145
146 for ( i = 0; i < N; i++ )
147 {
148 x[i] = 1.0;
149 }
150
151 for ( i = 0; i < N; i++ )
152 {
153 rhs[i] = 0.0;
154 for ( j = 0; j < N; j++ )
155 {
156 rhs[i] = rhs[i] + a[i+j*LDA] * x[j];
157 }
158 }
159
160 for ( i = 0; i < N; i++ )
161 {
162 resid[i] = -rhs[i];
163 for ( j = 0; j < N; j++ )
164 {
165 resid[i] = resid[i] + a[i+j*LDA] * b[j];
166 }
167 }
168
169 resid_max = 0.0;
170 for ( i = 0; i < N; i++ )
171 {
172 resid_max = r8_max ( resid_max, r8_abs ( resid[i] ) );
173 }
174
175 b_max = 0.0;
176 for ( i = 0; i < N; i++ )
177 {
178 b_max = r8_max ( b_max, r8_abs ( b[i] ) );
179 }
180
181 eps = r8_epsilon ( );
182
183 residn = resid_max / ( double ) N / a_max / b_max / eps;
184
185 time[2] = total;
186 if ( 0.0 < total )
187 {
188 time[3] = ops / ( 1.0E+06 * total );
189 }
190 else
191 {
192 time[3] = -1.0;
193 }
194 time[4] = 2.0 / time[3];
195 time[5] = total / cray;
196
197 printf ( "\n" );
198 printf ( " Norm. Resid Resid MACHEP X[1] X[N]\n" );
199 printf ( "\n" );
200 printf ( " %14f %14f %14e %14f %14f\n", residn, resid_max, eps, b[0], b[N-1] );
201 printf ( "\n" );
202 printf ( " Factor Solve Total Unit Cray-Ratio\n" );
203 printf ( "\n" );
204 printf ( " %9f %9f %9f %9f %9f\n",
205 time[0], time[1], time[2], time[4], time[5] );
206 printf ( "\n" );
207 printf ( "Unrolled Double Precision %9f Mflops\n", time[3]);
208 printf ( "\n" );
209
210 free ( a );
211 free ( b );
212 free ( ipvt );
213 free ( resid );
214 free ( rhs );
215 free ( x );
216 /*
217 Terminate.
218 */
219 printf ( "\n" );
220 printf ( "LINPACK_BENCH\n" );
221 printf ( " Normal end of execution.\n" );
222
223 printf ( "\n" );
224 timestamp ( );
225
226 return 0;
227 # undef LDA
228 # undef N
229 }
230 /******************************************************************************/
231
cpu_time(void)232 double cpu_time ( void )
233
234 /******************************************************************************/
235 /*
236 Purpose:
237
238 CPU_TIME returns the current reading on the CPU clock.
239
240 Discussion:
241
242 The CPU time measurements available through this routine are often
243 not very accurate. In some cases, the accuracy is no better than
244 a hundredth of a second.
245
246 Licensing:
247
248 This code is distributed under the GNU LGPL license.
249
250 Modified:
251
252 06 June 2005
253
254 Author:
255
256 John Burkardt
257
258 Parameters:
259
260 Output, double CPU_TIME, the current reading of the CPU clock, in seconds.
261 */
262 {
263 double value;
264
265 value = ( double ) clock ( )
266 / ( double ) CLOCKS_PER_SEC;
267
268 return value;
269 }
270 /******************************************************************************/
271
daxpy(int n,double da,double dx[],int incx,double dy[],int incy)272 void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy )
273
274 /******************************************************************************/
275 /*
276 Purpose:
277
278 DAXPY computes constant times a vector plus a vector.
279
280 Discussion:
281
282 This routine uses unrolled loops for increments equal to one.
283
284 Modified:
285
286 30 March 2007
287
288 Author:
289
290 FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
291 C version by John Burkardt
292
293 Reference:
294
295 Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
296 LINPACK User's Guide,
297 SIAM, 1979.
298
299 Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
300 Basic Linear Algebra Subprograms for Fortran Usage,
301 Algorithm 539,
302 ACM Transactions on Mathematical Software,
303 Volume 5, Number 3, September 1979, pages 308-323.
304
305 Parameters:
306
307 Input, int N, the number of elements in DX and DY.
308
309 Input, double DA, the multiplier of DX.
310
311 Input, double DX[*], the first vector.
312
313 Input, int INCX, the increment between successive entries of DX.
314
315 Input/output, double DY[*], the second vector.
316 On output, DY[*] has been replaced by DY[*] + DA * DX[*].
317
318 Input, int INCY, the increment between successive entries of DY.
319 */
320 {
321 int i;
322 int ix;
323 int iy;
324 int m;
325
326 if ( n <= 0 )
327 {
328 return;
329 }
330
331 if ( da == 0.0 )
332 {
333 return;
334 }
335 /*
336 Code for unequal increments or equal increments
337 not equal to 1.
338 */
339 if ( incx != 1 || incy != 1 )
340 {
341 if ( 0 <= incx )
342 {
343 ix = 0;
344 }
345 else
346 {
347 ix = ( - n + 1 ) * incx;
348 }
349
350 if ( 0 <= incy )
351 {
352 iy = 0;
353 }
354 else
355 {
356 iy = ( - n + 1 ) * incy;
357 }
358
359 for ( i = 0; i < n; i++ )
360 {
361 dy[iy] = dy[iy] + da * dx[ix];
362 ix = ix + incx;
363 iy = iy + incy;
364 }
365 }
366 /*
367 Code for both increments equal to 1.
368 */
369 else
370 {
371 m = n % 4;
372
373 for ( i = 0; i < m; i++ )
374 {
375 dy[i] = dy[i] + da * dx[i];
376 }
377
378 for ( i = m; i < n; i = i + 4 )
379 {
380 dy[i ] = dy[i ] + da * dx[i ];
381 dy[i+1] = dy[i+1] + da * dx[i+1];
382 dy[i+2] = dy[i+2] + da * dx[i+2];
383 dy[i+3] = dy[i+3] + da * dx[i+3];
384 }
385 }
386 return;
387 }
388 /******************************************************************************/
389
ddot(int n,double dx[],int incx,double dy[],int incy)390 double ddot ( int n, double dx[], int incx, double dy[], int incy )
391
392 /******************************************************************************/
393 /*
394 Purpose:
395
396 DDOT forms the dot product of two vectors.
397
398 Discussion:
399
400 This routine uses unrolled loops for increments equal to one.
401
402 Modified:
403
404 30 March 2007
405
406 Author:
407
408 FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
409 C version by John Burkardt
410
411 Reference:
412
413 Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
414 LINPACK User's Guide,
415 SIAM, 1979.
416
417 Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
418 Basic Linear Algebra Subprograms for Fortran Usage,
419 Algorithm 539,
420 ACM Transactions on Mathematical Software,
421 Volume 5, Number 3, September 1979, pages 308-323.
422
423 Parameters:
424
425 Input, int N, the number of entries in the vectors.
426
427 Input, double DX[*], the first vector.
428
429 Input, int INCX, the increment between successive entries in DX.
430
431 Input, double DY[*], the second vector.
432
433 Input, int INCY, the increment between successive entries in DY.
434
435 Output, double DDOT, the sum of the product of the corresponding
436 entries of DX and DY.
437 */
438 {
439 double dtemp;
440 int i;
441 int ix;
442 int iy;
443 int m;
444
445 dtemp = 0.0;
446
447 if ( n <= 0 )
448 {
449 return dtemp;
450 }
451 /*
452 Code for unequal increments or equal increments
453 not equal to 1.
454 */
455 if ( incx != 1 || incy != 1 )
456 {
457 if ( 0 <= incx )
458 {
459 ix = 0;
460 }
461 else
462 {
463 ix = ( - n + 1 ) * incx;
464 }
465
466 if ( 0 <= incy )
467 {
468 iy = 0;
469 }
470 else
471 {
472 iy = ( - n + 1 ) * incy;
473 }
474
475 for ( i = 0; i < n; i++ )
476 {
477 dtemp = dtemp + dx[ix] * dy[iy];
478 ix = ix + incx;
479 iy = iy + incy;
480 }
481 }
482 /*
483 Code for both increments equal to 1.
484 */
485 else
486 {
487 m = n % 5;
488
489 for ( i = 0; i < m; i++ )
490 {
491 dtemp = dtemp + dx[i] * dy[i];
492 }
493
494 for ( i = m; i < n; i = i + 5 )
495 {
496 dtemp = dtemp + dx[i ] * dy[i ]
497 + dx[i+1] * dy[i+1]
498 + dx[i+2] * dy[i+2]
499 + dx[i+3] * dy[i+3]
500 + dx[i+4] * dy[i+4];
501 }
502 }
503 return dtemp;
504 }
505 /******************************************************************************/
506
dgefa(double a[],int lda,int n,int ipvt[])507 int dgefa ( double a[], int lda, int n, int ipvt[] )
508
509 /******************************************************************************/
510 /*
511 Purpose:
512
513 DGEFA factors a real general matrix.
514
515 Modified:
516
517 16 May 2005
518
519 Author:
520
521 C version by John Burkardt.
522
523 Reference:
524
525 Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
526 LINPACK User's Guide,
527 SIAM, (Society for Industrial and Applied Mathematics),
528 3600 University City Science Center,
529 Philadelphia, PA, 19104-2688.
530 ISBN 0-89871-172-X
531
532 Parameters:
533
534 Input/output, double A[LDA*N].
535 On intput, the matrix to be factored.
536 On output, an upper triangular matrix and the multipliers used to obtain
537 it. The factorization can be written A=L*U, where L is a product of
538 permutation and unit lower triangular matrices, and U is upper triangular.
539
540 Input, int LDA, the leading dimension of A.
541
542 Input, int N, the order of the matrix A.
543
544 Output, int IPVT[N], the pivot indices.
545
546 Output, int DGEFA, singularity indicator.
547 0, normal value.
548 K, if U(K,K) == 0. This is not an error condition for this subroutine,
549 but it does indicate that DGESL or DGEDI will divide by zero if called.
550 Use RCOND in DGECO for a reliable indication of singularity.
551 */
552 {
553 int info;
554 int j;
555 int k;
556 int l;
557 double t;
558 /*
559 Gaussian elimination with partial pivoting.
560 */
561 info = 0;
562
563 for ( k = 1; k <= n-1; k++ )
564 {
565 /*
566 Find L = pivot index.
567 */
568 l = idamax ( n-k+1, a+(k-1)+(k-1)*lda, 1 ) + k - 1;
569 ipvt[k-1] = l;
570 /*
571 Zero pivot implies this column already triangularized.
572 */
573 if ( a[l-1+(k-1)*lda] == 0.0 )
574 {
575 info = k;
576 continue;
577 }
578 /*
579 Interchange if necessary.
580 */
581 if ( l != k )
582 {
583 t = a[l-1+(k-1)*lda];
584 a[l-1+(k-1)*lda] = a[k-1+(k-1)*lda];
585 a[k-1+(k-1)*lda] = t;
586 }
587 /*
588 Compute multipliers.
589 */
590 t = -1.0 / a[k-1+(k-1)*lda];
591
592 dscal ( n-k, t, a+k+(k-1)*lda, 1 );
593 /*
594 Row elimination with column indexing.
595 */
596 for ( j = k+1; j <= n; j++ )
597 {
598 t = a[l-1+(j-1)*lda];
599 if ( l != k )
600 {
601 a[l-1+(j-1)*lda] = a[k-1+(j-1)*lda];
602 a[k-1+(j-1)*lda] = t;
603 }
604 daxpy ( n-k, t, a+k+(k-1)*lda, 1, a+k+(j-1)*lda, 1 );
605 }
606
607 }
608
609 ipvt[n-1] = n;
610
611 if ( a[n-1+(n-1)*lda] == 0.0 )
612 {
613 info = n;
614 }
615
616 return info;
617 }
618 /******************************************************************************/
619
dgesl(double a[],int lda,int n,int ipvt[],double b[],int job)620 void dgesl ( double a[], int lda, int n, int ipvt[], double b[], int job )
621
622 /******************************************************************************/
623 /*
624 Purpose:
625
626 DGESL solves a real general linear system A * X = B.
627
628 Discussion:
629
630 DGESL can solve either of the systems A * X = B or A' * X = B.
631
632 The system matrix must have been factored by DGECO or DGEFA.
633
634 A division by zero will occur if the input factor contains a
635 zero on the diagonal. Technically this indicates singularity
636 but it is often caused by improper arguments or improper
637 setting of LDA. It will not occur if the subroutines are
638 called correctly and if DGECO has set 0.0 < RCOND
639 or DGEFA has set INFO == 0.
640
641 Modified:
642
643 16 May 2005
644
645 Author:
646
647 C version by John Burkardt.
648
649 Reference:
650
651 Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
652 LINPACK User's Guide,
653 SIAM, (Society for Industrial and Applied Mathematics),
654 3600 University City Science Center,
655 Philadelphia, PA, 19104-2688.
656 ISBN 0-89871-172-X
657
658 Parameters:
659
660 Input, double A[LDA*N], the output from DGECO or DGEFA.
661
662 Input, int LDA, the leading dimension of A.
663
664 Input, int N, the order of the matrix A.
665
666 Input, int IPVT[N], the pivot vector from DGECO or DGEFA.
667
668 Input/output, double B[N].
669 On input, the right hand side vector.
670 On output, the solution vector.
671
672 Input, int JOB.
673 0, solve A * X = B;
674 nonzero, solve A' * X = B.
675 */
676 {
677 int k;
678 int l;
679 double t;
680 /*
681 Solve A * X = B.
682 */
683 if ( job == 0 )
684 {
685 for ( k = 1; k <= n-1; k++ )
686 {
687 l = ipvt[k-1];
688 t = b[l-1];
689
690 if ( l != k )
691 {
692 b[l-1] = b[k-1];
693 b[k-1] = t;
694 }
695
696 daxpy ( n-k, t, a+k+(k-1)*lda, 1, b+k, 1 );
697
698 }
699
700 for ( k = n; 1 <= k; k-- )
701 {
702 b[k-1] = b[k-1] / a[k-1+(k-1)*lda];
703 t = -b[k-1];
704 daxpy ( k-1, t, a+0+(k-1)*lda, 1, b, 1 );
705 }
706 }
707 /*
708 Solve A' * X = B.
709 */
710 else
711 {
712 for ( k = 1; k <= n; k++ )
713 {
714 t = ddot ( k-1, a+0+(k-1)*lda, 1, b, 1 );
715 b[k-1] = ( b[k-1] - t ) / a[k-1+(k-1)*lda];
716 }
717
718 for ( k = n-1; 1 <= k; k-- )
719 {
720 b[k-1] = b[k-1] + ddot ( n-k, a+k+(k-1)*lda, 1, b+k, 1 );
721 l = ipvt[k-1];
722
723 if ( l != k )
724 {
725 t = b[l-1];
726 b[l-1] = b[k-1];
727 b[k-1] = t;
728 }
729 }
730 }
731 return;
732 }
733 /******************************************************************************/
734
dscal(int n,double sa,double x[],int incx)735 void dscal ( int n, double sa, double x[], int incx )
736
737 /******************************************************************************/
738 /*
739 Purpose:
740
741 DSCAL scales a vector by a constant.
742
743 Modified:
744
745 30 March 2007
746
747 Author:
748
749 FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
750 C version by John Burkardt
751
752 Reference:
753
754 Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
755 LINPACK User's Guide,
756 SIAM, 1979.
757
758 Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
759 Basic Linear Algebra Subprograms for Fortran Usage,
760 Algorithm 539,
761 ACM Transactions on Mathematical Software,
762 Volume 5, Number 3, September 1979, pages 308-323.
763
764 Parameters:
765
766 Input, int N, the number of entries in the vector.
767
768 Input, double SA, the multiplier.
769
770 Input/output, double X[*], the vector to be scaled.
771
772 Input, int INCX, the increment between successive entries of X.
773 */
774 {
775 int i;
776 int ix;
777 int m;
778
779 if ( n <= 0 )
780 {
781 }
782 else if ( incx == 1 )
783 {
784 m = n % 5;
785
786 for ( i = 0; i < m; i++ )
787 {
788 x[i] = sa * x[i];
789 }
790
791 for ( i = m; i < n; i = i + 5 )
792 {
793 x[i] = sa * x[i];
794 x[i+1] = sa * x[i+1];
795 x[i+2] = sa * x[i+2];
796 x[i+3] = sa * x[i+3];
797 x[i+4] = sa * x[i+4];
798 }
799 }
800 else
801 {
802 if ( 0 <= incx )
803 {
804 ix = 0;
805 }
806 else
807 {
808 ix = ( - n + 1 ) * incx;
809 }
810
811 for ( i = 0; i < n; i++ )
812 {
813 x[ix] = sa * x[ix];
814 ix = ix + incx;
815 }
816 }
817 return;
818 }
819 /******************************************************************************/
820
idamax(int n,double dx[],int incx)821 int idamax ( int n, double dx[], int incx )
822
823 /******************************************************************************/
824 /*
825 Purpose:
826
827 IDAMAX finds the index of the vector element of maximum absolute value.
828
829 Discussion:
830
831 WARNING: This index is a 1-based index, not a 0-based index!
832
833 Modified:
834
835 30 March 2007
836
837 Author:
838
839 FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
840 C version by John Burkardt
841
842 Reference:
843
844 Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
845 LINPACK User's Guide,
846 SIAM, 1979.
847
848 Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
849 Basic Linear Algebra Subprograms for Fortran Usage,
850 Algorithm 539,
851 ACM Transactions on Mathematical Software,
852 Volume 5, Number 3, September 1979, pages 308-323.
853
854 Parameters:
855
856 Input, int N, the number of entries in the vector.
857
858 Input, double X[*], the vector to be examined.
859
860 Input, int INCX, the increment between successive entries of SX.
861
862 Output, int IDAMAX, the index of the element of maximum
863 absolute value.
864 */
865 {
866 double dmax;
867 int i;
868 int ix;
869 int value;
870
871 value = 0;
872
873 if ( n < 1 || incx <= 0 )
874 {
875 return value;
876 }
877
878 value = 1;
879
880 if ( n == 1 )
881 {
882 return value;
883 }
884
885 if ( incx == 1 )
886 {
887 dmax = r8_abs ( dx[0] );
888
889 for ( i = 1; i < n; i++ )
890 {
891 if ( dmax < r8_abs ( dx[i] ) )
892 {
893 value = i + 1;
894 dmax = r8_abs ( dx[i] );
895 }
896 }
897 }
898 else
899 {
900 ix = 0;
901 dmax = r8_abs ( dx[0] );
902 ix = ix + incx;
903
904 for ( i = 1; i < n; i++ )
905 {
906 if ( dmax < r8_abs ( dx[ix] ) )
907 {
908 value = i + 1;
909 dmax = r8_abs ( dx[ix] );
910 }
911 ix = ix + incx;
912 }
913 }
914
915 return value;
916 }
917 /******************************************************************************/
918
r8_abs(double x)919 double r8_abs ( double x )
920
921 /******************************************************************************/
922 /*
923 Purpose:
924
925 R8_ABS returns the absolute value of a R8.
926
927 Modified:
928
929 02 April 2005
930
931 Author:
932
933 John Burkardt
934
935 Parameters:
936
937 Input, double X, the quantity whose absolute value is desired.
938
939 Output, double R8_ABS, the absolute value of X.
940 */
941 {
942 double value;
943
944 if ( 0.0 <= x )
945 {
946 value = x;
947 }
948 else
949 {
950 value = -x;
951 }
952 return value;
953 }
954 /******************************************************************************/
955
r8_epsilon()956 double r8_epsilon ( )
957
958 /******************************************************************************/
959 /*
960 Purpose:
961
962 R8_EPSILON returns the R8 round off unit.
963
964 Discussion:
965
966 R8_EPSILON is a number R which is a power of 2 with the property that,
967 to the precision of the computer's arithmetic,
968 1 < 1 + R
969 but
970 1 = ( 1 + R / 2 )
971
972 Licensing:
973
974 This code is distributed under the GNU LGPL license.
975
976 Modified:
977
978 01 September 2012
979
980 Author:
981
982 John Burkardt
983
984 Parameters:
985
986 Output, double R8_EPSILON, the R8 round-off unit.
987 */
988 {
989 const double value = 2.220446049250313E-016;
990
991 return value;
992 }
993 /******************************************************************************/
994
r8_max(double x,double y)995 double r8_max ( double x, double y )
996
997 /******************************************************************************/
998 /*
999 Purpose:
1000
1001 R8_MAX returns the maximum of two R8's.
1002
1003 Modified:
1004
1005 18 August 2004
1006
1007 Author:
1008
1009 John Burkardt
1010
1011 Parameters:
1012
1013 Input, double X, Y, the quantities to compare.
1014
1015 Output, double R8_MAX, the maximum of X and Y.
1016 */
1017 {
1018 double value;
1019
1020 if ( y < x )
1021 {
1022 value = x;
1023 }
1024 else
1025 {
1026 value = y;
1027 }
1028 return value;
1029 }
1030 /******************************************************************************/
1031
r8_random(int iseed[4])1032 double r8_random ( int iseed[4] )
1033
1034 /******************************************************************************/
1035 /*
1036 Purpose:
1037
1038 R8_RANDOM returns a uniformly distributed random number between 0 and 1.
1039
1040 Discussion:
1041
1042 This routine uses a multiplicative congruential method with modulus
1043 2**48 and multiplier 33952834046453 (see G.S.Fishman,
1044 'Multiplicative congruential random number generators with modulus
1045 2**b: an exhaustive analysis for b = 32 and a partial analysis for
1046 b = 48', Math. Comp. 189, pp 331-344, 1990).
1047
1048 48-bit integers are stored in 4 integer array elements with 12 bits
1049 per element. Hence the routine is portable across machines with
1050 integers of 32 bits or more.
1051
1052 Parameters:
1053
1054 Input/output, integer ISEED(4).
1055 On entry, the seed of the random number generator; the array
1056 elements must be between 0 and 4095, and ISEED(4) must be odd.
1057 On exit, the seed is updated.
1058
1059 Output, double R8_RANDOM, the next pseudorandom number.
1060 */
1061 {
1062 int ipw2 = 4096;
1063 int it1;
1064 int it2;
1065 int it3;
1066 int it4;
1067 int m1 = 494;
1068 int m2 = 322;
1069 int m3 = 2508;
1070 int m4 = 2549;
1071 double one = 1.0;
1072 double r = 1.0 / 4096.0;
1073 double value;
1074 /*
1075 Multiply the seed by the multiplier modulo 2**48.
1076 */
1077 it4 = iseed[3] * m4;
1078 it3 = it4 / ipw2;
1079 it4 = it4 - ipw2 * it3;
1080 it3 = it3 + iseed[2] * m4 + iseed[3] * m3;
1081 it2 = it3 / ipw2;
1082 it3 = it3 - ipw2 * it2;
1083 it2 = it2 + iseed[1] * m4 + iseed[2] * m3 + iseed[3] * m2;
1084 it1 = it2 / ipw2;
1085 it2 = it2 - ipw2 * it1;
1086 it1 = it1 + iseed[0] * m4 + iseed[1] * m3 + iseed[2] * m2 + iseed[3] * m1;
1087 it1 = ( it1 % ipw2 );
1088 /*
1089 Return updated seed
1090 */
1091 iseed[0] = it1;
1092 iseed[1] = it2;
1093 iseed[2] = it3;
1094 iseed[3] = it4;
1095 /*
1096 Convert 48-bit integer to a real number in the interval (0,1)
1097 */
1098 value =
1099 r * ( ( double ) ( it1 )
1100 + r * ( ( double ) ( it2 )
1101 + r * ( ( double ) ( it3 )
1102 + r * ( ( double ) ( it4 ) ) ) ) );
1103
1104 return value;
1105 }
1106 /******************************************************************************/
1107
r8mat_gen(int lda,int n)1108 double *r8mat_gen ( int lda, int n )
1109
1110 /******************************************************************************/
1111 /*
1112 Purpose:
1113
1114 R8MAT_GEN generates a random R8MAT.
1115
1116 Modified:
1117
1118 06 June 2005
1119
1120 Parameters:
1121
1122 Input, integer LDA, the leading dimension of the matrix.
1123
1124 Input, integer N, the order of the matrix.
1125
1126 Output, double R8MAT_GEN[LDA*N], the N by N matrix.
1127 */
1128 {
1129 double *a;
1130 int i;
1131 int init[4] = { 1, 2, 3, 1325 };
1132 int j;
1133
1134 a = ( double * ) malloc ( lda * n * sizeof ( double ) );
1135
1136 for ( j = 1; j <= n; j++ )
1137 {
1138 for ( i = 1; i <= n; i++ )
1139 {
1140 a[i-1+(j-1)*lda] = r8_random ( init ) - 0.5;
1141 }
1142 }
1143
1144 return a;
1145 }
1146 /******************************************************************************/
1147
timestamp(void)1148 void timestamp ( void )
1149
1150 /******************************************************************************/
1151 /*
1152 Purpose:
1153
1154 TIMESTAMP prints the current YMDHMS date as a time stamp.
1155
1156 Example:
1157
1158 31 May 2001 09:45:54 AM
1159
1160 Licensing:
1161
1162 This code is distributed under the GNU LGPL license.
1163
1164 Modified:
1165
1166 24 September 2003
1167
1168 Author:
1169
1170 John Burkardt
1171
1172 Parameters:
1173
1174 None
1175 */
1176 {
1177 # define TIME_SIZE 40
1178
1179 static char time_buffer[TIME_SIZE];
1180 const struct tm *tm;
1181 size_t len;
1182 time_t now;
1183
1184 now = time ( NULL );
1185 tm = localtime ( &now );
1186
1187 len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
1188
1189 printf ( "%s\n", time_buffer );
1190
1191 return;
1192 # undef TIME_SIZE
1193 }
1194
1195