Jovan Ivković
/
Linpack
Important changes to repositories hosted on mbed.com
Mbed hosted mercurial repositories are deprecated and are due to be permanently deleted in July 2026.
To keep a copy of this software download the repository Zip archive or clone locally using Mercurial.
It is also possible to export all your personal repositories from the account settings page.
Diff: main.cpp
- Revision:
- 2:273153e44338
- Parent:
- 1:be78b18b8347
- Child:
- 3:da1132c65314
--- a/main.cpp Mon Jan 02 02:55:50 2017 +0000
+++ b/main.cpp Tue Jan 03 02:31:38 2017 +0000
@@ -1,424 +1,1108 @@
//********************************************************
-//** Nucleo-144 Stm32F746 and Stm32F767 benchmark ******
+//** Nucleo-144 Stm32F746 and Stm32F767 benchmark ******
+//** Limpack -port form Arduino IDE *****
//** Jovan Ivkovic - 2016 ******
//********************************************************
#include "mbed.h"
+
+/* the following is optional depending on the timing function used */
+# include <stdlib.h>
+# include <stdio.h>
+# include <math.h>
+
DigitalOut myled(LED1);
Serial pc(USBTX, USBRX);
Timer timer;
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <math.h>
-/* the following is optional depending on the timing function used */
-#include <time.h>
+int do_benchmark( void );
+double cpu_time( void );
+void daxpy( int n, double da, double dx[], int incx, double dy[], int incy );
+double ddot( int n, double dx[], int incx, double dy[], int incy );
+int dgefa( double a[], int lda, int n, int ipvt[] );
+void dgesl( double a[], int lda, int n, int ipvt[], double b[], int job );
+void dscal( int n, double sa, double x[], int incx );
+int idamax( int n, double dx[], int incx );
+double r8_abs( double x );
+double r8_epsilon( void );
+double r8_max( double x, double y );
+double r8_random(int iseed[4] );
+double *r8mat_gen ( int lda, int n );
-/* map the FORTRAN math functions, etc. to the C versions */
-#define DSIN sin
-#define DCOS cos
-#define DATAN atan
-#define DLOG log
-#define DEXP exp
-#define DSQRT sqrt
-#define IF if
+//static FILE uartout = {0} ;
-/* function prototypes */
-void POUT(long N, long J, long K, double X1, double X2, double X3, double X4);
-void PA(double E[]);
-void P0(void);
-void P3(double X, double Y, double *Z);
-#define USAGE "usage: whetdc [-c] [loops]\n"
+//static int uart_putchar (char c, FILE *stream)
+//{
+// Serial.write(c) ;
+// return 0 ;
+//}
-/*
- COMMON T,T1,T2,E1(4),J,K,L
-*/
-double T,T1,T2,E1[5];
-int J,K,L;
-int argc = 0; //Mod for nucleo. Change in code below if you want non-default loop count
+//void setup() {
+// Serial.begin(9600);
+// fdev_setup_stream (&uartout, uart_putchar, NULL, _FDEV_SETUP_WRITE);
+// stdout = &uartout ;
+//}
-//************************************
-//** Whetstone 64b-DP **
-//** SUB **
-//************************************
-int Whetstone() // ------------ Metoda -----------
+int main()
{
- pc.baud(115200);
- pc.printf("Beginning Whetstone benchmark at ");
-
- pc.printf("default 216 MHz ...\n");
- /* used in the FORTRAN version */
- long I;
- long N1, N2, N3, N4, N6, N7, N8, N9, N10, N11;
- double X1,X2,X3,X4,X,Y,Z;
- long LOOP;
- int II, JJ;
+ //pc.baud(115200);
+ pc.baud(19200);
+
+ while(1) {
+
+ pc.printf("Starting benchmark...\n");
+
+ do_benchmark();
+
+ pc.printf(" kraj \n\n");
+ }
+}
+
+/******************************************************************************/
- /* added for this version */
- long loopstart = 0;
- long startsec,finisec = 0;
- double KIPS;
- int continuous;
+int do_benchmark ( void )
- loopstart = 25000; /* 1000 see the note about LOOP below */
- continuous = 0;
+/******************************************************************************/
+/*
+ Purpose:
+
+ MAIN is the main program for LINPACK_BENCH.
- II = 1; /* start at the first arg (temp use of II here) */
-
-LCONT:
-/*
-********************************************
-* Start benchmark timing at this point.
-********************************************
-*/
- timer.start();
- startsec = 0;
- finisec = 0;
- startsec = timer.read_us();
+ Discussion:
+
+ LINPACK_BENCH drives the double precision LINPACK benchmark program.
+
+ Modified:
+
+ 25 July 2008
+
+ Parameters:
-/*
-********************************************
-* The actual benchmark starts here.
-********************************************
+ N is the problem size.
*/
- T = .499975;
- T1 = 0.50025;
- T2 = 2.0;
-/*
-********************************************
-* With loopcount LOOP=10, one million Whetstone instructions
-* will be executed in EACH MAJOR LOOP..A MAJOR LOOP IS EXECUTED
-* 'II' TIMES TO INCREASE WALL-CLOCK TIMING ACCURACY.
-*
-* LOOP = 1000;
-*/
- LOOP = loopstart;
- II = 1;
- JJ = 1;
+{
+# define N 8
+# define LDA ( N + 1 )
+
+ static double *a;
+ static double a_max;
+ static double *b;
+ static double b_max;
+ const double cray = 0.056;
+ static double eps;
+ int i;
+ int info;
+ static int ipvt[N];
+ int j;
+ int job;
+ double ops;
+ static double resid[N];
+ double resid_max;
+ double residn;
+ static double rhs[N];
+ float t1;
+ float t2;
+ static double time[6];
+ double total;
+ double x[N];
+
-IILOOP:
- N1 = 0;
- N2 = 12 * LOOP;
- N3 = 14 * LOOP;
- N4 = 345 * LOOP;
- N6 = 210 * LOOP;
- N7 = 32 * LOOP;
- N8 = 899 * LOOP;
- N9 = 616 * LOOP;
- N10 = 0;
- N11 = 93 * LOOP;
-/*
-********************************************
-* Module 1: Simple identifiers
-********************************************
-*/
- X1 = 1.0;
- X2 = -1.0;
- X3 = -1.0;
- X4 = -1.0;
+ pc.printf ( "\n" );
+ pc.printf ( "LINPACK_BENCH\n" );
+ pc.printf ( " C version\n" );
+ pc.printf ( "\n" );
+ pc.printf ( " The LINPACK benchmark.\n" );
+ pc.printf ( " Language: C\n" );
+ pc.printf ( " Datatype: Double precision real\n" );
+ pc.printf ( " Matrix order N = %d\n", N );
+ pc.printf ( " Leading matrix dimension LDA = %d\n", LDA );
+
+ ops = ( double ) ( 2L * N * N * N ) / 3.0 + 2.0 * ( double ) ( ( long ) N * N );
+
+ /*
+ Allocate space for arrays.
+ */
+ a = r8mat_gen ( LDA, N );
- for (I = 1; I <= N1; I++)
- {
- X1 = (X1 + X2 + X3 - X4) * T;
- X2 = (X1 + X2 - X3 + X4) * T;
- X3 = (X1 - X2 + X3 + X4) * T;
- X4 = (-X1+ X2 + X3 + X4) * T;
+ a_max = 0.0;
+ for ( j = 0; j < N; j++ ) {
+ for ( i = 0; i < N; i++ ) {
+ a_max = r8_max ( a_max, a[i+j*LDA] );
+ }
}
-#ifdef PRINTOUT
- IF (JJ==II) POUT(N1,N1,N1,X1,X2,X3,X4);
-#endif
-/*
-********************************************
-* Module 2: Array elements
-********************************************
-*/
- E1[1] = 1.0;
- E1[2] = -1.0;
- E1[3] = -1.0;
- E1[4] = -1.0;
+ for ( i = 0; i < N; i++ ) {
+ x[i] = 1.0;
+ }
- for (I = 1; I <= N2; I++)
- {
- E1[1] = ( E1[1] + E1[2] + E1[3] - E1[4]) * T;
- E1[2] = ( E1[1] + E1[2] - E1[3] + E1[4]) * T;
- E1[3] = ( E1[1] - E1[2] + E1[3] + E1[4]) * T;
- E1[4] = (-E1[1] + E1[2] + E1[3] + E1[4]) * T;
+ for ( i = 0; i < N; i++ ) {
+ b[i] = 0.0;
+ for ( j = 0; j < N; j++ ) {
+ b[i] = b[i] + a[i+j*LDA] * x[j];
+ }
}
-#ifdef PRINTOUT
- IF (JJ==II) POUT(N2,N3,N2,E1[1],E1[2],E1[3],E1[4]);
-#endif
+ timer.start();
+
+ //*****************
+ t1 = timer.read();
+
+ info = dgefa ( a, LDA, N, ipvt );
-/*
-********************************************
-* Module 3: Array as parameter
-********************************************
-*/
- for (I = 1; I <= N3; I++)
- {
- PA(E1);
+ t2 = timer.read();
+
+ if ( info != 0 ) {
+ pc.printf ( "\n" );
+ pc.printf ( "LINPACK_BENCH - Fatal error!\n" );
+ pc.printf ( " The matrix A is apparently singular.\n" );
+ pc.printf ( " Abnormal end of execution.\n" );
+ return 1;
}
-#ifdef PRINTOUT
- IF (JJ==II) POUT(N3,N2,N2,E1[1],E1[2],E1[3],E1[4]);
-#endif
+ time[0] = ( double )(t2 - t1);
+
+ //*********
+ timer.reset();
+ t1 = timer.read;
+
+ job = 0;
+ dgesl ( a, LDA, N, ipvt, b, job );
+
+ t2 = timer.read;
+ time[1] = ( double )(t2 - t1);
+
+ total = time[0] + time[1];
-/*
-********************************************
-* Module 4: Conditional jumps
-********************************************
-*/
- J = 1;
- for (I = 1; I <= N4; I++)
- {
- if (J == 1)
- J = 2;
- else
- J = 3;
+ timer.reset();
+ //*********
+
+ /*
+ Compute a residual to verify results.
+ */
+ a = r8mat_gen ( LDA, N );
+
+ for ( i = 0; i < N; i++ ) {
+ x[i] = 1.0;
+ }
- if (J > 2)
- J = 0;
- else
- J = 1;
+ for ( i = 0; i < N; i++ ) {
+ rhs[i] = 0.0;
+ for ( j = 0; j < N; j++ ) {
+ rhs[i] = rhs[i] + a[i+j*LDA] * x[j];
+ }
+ }
- if (J < 1)
- J = 1;
- else
- J = 0;
+ for ( i = 0; i < N; i++ ) {
+ resid[i] = -rhs[i];
+ for ( j = 0; j < N; j++ ) {
+ resid[i] = resid[i] + a[i+j*LDA] * b[j];
+ }
+ }
+
+ resid_max = 0.0;
+ for ( i = 0; i < N; i++ ) {
+ resid_max = r8_max ( resid_max, r8_abs ( resid[i] ) );
+ }
+
+ b_max = 0.0;
+ for ( i = 0; i < N; i++ ) {
+ b_max = r8_max ( b_max, r8_abs ( b[i] ) );
}
-#ifdef PRINTOUT
- IF (JJ==II) POUT(N4,J,J,X1,X2,X3,X4);
-#endif
+ eps = r8_epsilon ( );
+
+ residn = resid_max / ( double ) N / a_max / b_max / eps;
+
+ time[2] = total;
+
+ if ( total > 0.0 )
+ {
+ time[3] = ops / ( 1.0E+06 * total );
+ }
+ else
+ {
+ time[3] = -1.0;
+ }
+
+ time[4] = 2.0 / time[3];
+ time[5] = total / cray;
+
+ pc.printf( "\n" );
+ pc.printf( " Norm. Resid Resid MACHEP X[1] X[N] \n\n" );
+ /* pc.printf(" %14f", residn);
+
+ pc.printf(" %14f", resid_max);
+
+ pc.printf(" %f14", eps);
+
+ pc.printf(" %f14", b[0]);
+ pc.printf(" %14f ",b[N-1]);
+ pc.printf("\n\n");
+ */
+ //pc.printf( " %14f %14f %14e %14f %14f \n", residn, resid_max, eps, b[0], b[N-1] );
+
+ pc.printf( " \n\n ");
+
+ pc.printf( " Factor Solve Total MFLOPS Unit Cray-Ratio \n\n" );
+ // */
+ for(int ii=0; ii<6; ii++) {
+ pc.printf("\t %9f", time[ii]);
+ }
+
+ // pc.printf( " %9f %9f %9f %9f %9f %9f\n", time[0], time[1], time[2], time[3], time[4], time[5] );
+
+ /*
+ Terminate.
+ */
+ pc.printf( "\n" );
+ pc.printf( "LINPACK_BENCH\n" );
+ pc.printf( " Normal end of execution.\n" );
+
+ pc.printf( "\n" );
+
+ timer.stop();
+
+ return 0;
+# undef LDA
+# undef N
+}
+/******************************************************************************/
+
+double cpu_time ( void )
+
+/******************************************************************************/
/*
-********************************************
-* Module 5: Omitted
-* Module 6: Integer arithmetic
-********************************************
-*/
+ Purpose:
+
+ CPU_TIME returns the current reading on the CPU clock.
+
+ Discussion:
- J = 1;
- K = 2;
- L = 3;
+ The CPU time measurements available through this routine are often
+ not very accurate. In some cases, the accuracy is no better than
+ a hundredth of a second.
+
+ koristi mbed.Timer
+
+*/
+{
+ double vreme;
+
+ vreme = timer.read_ms() / 1000;
+
+ return vreme;
+}
+/******************************************************************************/
+
+void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ DAXPY computes constant times a vector plus a vector.
+
+ Discussion:
+
+ This routine uses unrolled loops for increments equal to one.
+
+ Modified:
+
+ 30 March 2007
+
+ Author:
- for (I = 1; I <= N6; I++)
- {
- J = J * (K-J) * (L-K);
- K = L * K - (L-J) * K;
- L = (L-K) * (K+J);
- E1[L-1] = J + K + L;
- E1[K-1] = J * K * L;
+ FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
+ C version by John Burkardt
+
+ Reference:
+
+ Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
+ LINPACK User's Guide,
+ SIAM, 1979.
+
+ Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
+ Basic Linear Algebra Subprograms for Fortran Usage,
+ Algorithm 539,
+ ACM Transactions on Mathematical Software,
+ Volume 5, Number 3, September 1979, pages 308-323.
+
+ Parameters:
+
+ Input, int N, the number of elements in DX and DY.
+
+ Input, double DA, the multiplier of DX.
+
+ Input, double DX[*], the first vector.
+
+ Input, int INCX, the increment between successive entries of DX.
+
+ Input/output, double DY[*], the second vector.
+ On output, DY[*] has been replaced by DY[*] + DA * DX[*].
+
+ Input, int INCY, the increment between successive entries of DY.
+*/
+{
+ int i;
+ int ix;
+ int iy;
+ int m;
+
+ if ( n <= 0 ) {
+ return;
}
-#ifdef PRINTOUT
- IF (JJ==II) POUT(N6,J,K,E1[1],E1[2],E1[3],E1[4]);
-#endif
+ if ( da == 0.0 ) {
+ return;
+ }
+ /*
+ Code for unequal increments or equal increments
+ not equal to 1.
+ */
+ if ( incx != 1 || incy != 1 ) {
+ if ( 0 <= incx ) {
+ ix = 0;
+ } else {
+ ix = ( - n + 1 ) * incx;
+ }
+
+ if ( 0 <= incy ) {
+ iy = 0;
+ } else {
+ iy = ( - n + 1 ) * incy;
+ }
+ for ( i = 0; i < n; i++ ) {
+ dy[iy] = dy[iy] + da * dx[ix];
+ ix = ix + incx;
+ iy = iy + incy;
+ }
+ }
+ /*
+ Code for both increments equal to 1.
+ */
+ else {
+ m = n % 4;
+
+ for ( i = 0; i < m; i++ ) {
+ dy[i] = dy[i] + da * dx[i];
+ }
+
+ for ( i = m; i < n; i = i + 4 ) {
+ dy[i ] = dy[i ] + da * dx[i ];
+ dy[i+1] = dy[i+1] + da * dx[i+1];
+ dy[i+2] = dy[i+2] + da * dx[i+2];
+ dy[i+3] = dy[i+3] + da * dx[i+3];
+ }
+ }
+ return;
+}
+/******************************************************************************/
+
+double ddot ( int n, double dx[], int incx, double dy[], int incy )
+
+/******************************************************************************/
/*
-********************************************
-* Module 7: Trigonometric functions
-********************************************
+ Purpose:
+
+ DDOT forms the dot product of two vectors.
+
+ Discussion:
+
+ This routine uses unrolled loops for increments equal to one.
+
+ Modified:
+
+ 30 March 2007
+
+ Author:
+
+ FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
+ C version by John Burkardt
+
+ Reference:
+
+ Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
+ LINPACK User's Guide,
+ SIAM, 1979.
+
+ Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
+ Basic Linear Algebra Subprograms for Fortran Usage,
+ Algorithm 539,
+ ACM Transactions on Mathematical Software,
+ Volume 5, Number 3, September 1979, pages 308-323.
+
+ Parameters:
+
+ Input, int N, the number of entries in the vectors.
+
+ Input, double DX[*], the first vector.
+
+ Input, int INCX, the increment between successive entries in DX.
+
+ Input, double DY[*], the second vector.
+
+ Input, int INCY, the increment between successive entries in DY.
+
+ Output, double DDOT, the sum of the product of the corresponding
+ entries of DX and DY.
*/
- X = 0.5;
- Y = 0.5;
+{
+ double dtemp;
+ int i;
+ int ix;
+ int iy;
+ int m;
+
+ dtemp = 0.0;
+
+ if ( n <= 0 ) {
+ return dtemp;
+ }
+ /*
+ Code for unequal increments or equal increments
+ not equal to 1.
+ */
+ if ( incx != 1 || incy != 1 ) {
+ if ( 0 <= incx ) {
+ ix = 0;
+ } else {
+ ix = ( - n + 1 ) * incx;
+ }
+
+ if ( 0 <= incy ) {
+ iy = 0;
+ } else {
+ iy = ( - n + 1 ) * incy;
+ }
+
+ for ( i = 0; i < n; i++ ) {
+ dtemp = dtemp + dx[ix] * dy[iy];
+ ix = ix + incx;
+ iy = iy + incy;
+ }
+ }
+ /*
+ Code for both increments equal to 1.
+ */
+ else {
+ m = n % 5;
+
+ for ( i = 0; i < m; i++ ) {
+ dtemp = dtemp + dx[i] * dy[i];
+ }
+
+ for ( i = m; i < n; i = i + 5 ) {
+ 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];
+ }
+ }
+ return dtemp;
+}
+/******************************************************************************/
+
+int dgefa ( double a[], int lda, int n, int ipvt[] )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ DGEFA factors a real general matrix.
+
+ Modified:
+
+ 16 May 2005
+
+ Author:
+
+ C version by John Burkardt.
+
+ Reference:
- for (I = 1; I <= N7; I++)
- {
- X = T * DATAN(T2*DSIN(X)*DCOS(X)/(DCOS(X+Y)+DCOS(X-Y)-1.0));
- Y = T * DATAN(T2*DSIN(Y)*DCOS(Y)/(DCOS(X+Y)+DCOS(X-Y)-1.0));
+ Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
+ LINPACK User's Guide,
+ SIAM, (Society for Industrial and Applied Mathematics),
+ 3600 University City Science Center,
+ Philadelphia, PA, 19104-2688.
+ ISBN 0-89871-172-X
+
+ Parameters:
+
+ Input/output, double A[LDA*N].
+ On intput, the matrix to be factored.
+ On output, an upper triangular matrix and the multipliers used to obtain
+ it. The factorization can be written A=L*U, where L is a product of
+ permutation and unit lower triangular matrices, and U is upper triangular.
+
+ Input, int LDA, the leading dimension of A.
+
+ Input, int N, the order of the matrix A.
+
+ Output, int IPVT[N], the pivot indices.
+
+ Output, int DGEFA, singularity indicator.
+ 0, normal value.
+ K, if U(K,K) == 0. This is not an error condition for this subroutine,
+ but it does indicate that DGESL or DGEDI will divide by zero if called.
+ Use RCOND in DGECO for a reliable indication of singularity.
+*/
+{
+ int info;
+ int j;
+ int k;
+ int l;
+ double t;
+ /*
+ Gaussian elimination with partial pivoting.
+ */
+ info = 0;
+
+ for ( k = 1; k <= n-1; k++ ) {
+ /*
+ Find L = pivot index.
+ */
+ l = idamax ( n-k+1, a+(k-1)+(k-1)*lda, 1 ) + k - 1;
+ ipvt[k-1] = l;
+ /*
+ Zero pivot implies this column already triangularized.
+ */
+ if ( a[l-1+(k-1)*lda] == 0.0 ) {
+ info = k;
+ continue;
+ }
+ /*
+ Interchange if necessary.
+ */
+ if ( l != k ) {
+ t = a[l-1+(k-1)*lda];
+ a[l-1+(k-1)*lda] = a[k-1+(k-1)*lda];
+ a[k-1+(k-1)*lda] = t;
+ }
+ /*
+ Compute multipliers.
+ */
+ t = -1.0 / a[k-1+(k-1)*lda];
+
+ dscal ( n-k, t, a+k+(k-1)*lda, 1 );
+ /*
+ Row elimination with column indexing.
+ */
+ for ( j = k+1; j <= n; j++ ) {
+ t = a[l-1+(j-1)*lda];
+ if ( l != k ) {
+ a[l-1+(j-1)*lda] = a[k-1+(j-1)*lda];
+ a[k-1+(j-1)*lda] = t;
+ }
+ daxpy ( n-k, t, a+k+(k-1)*lda, 1, a+k+(j-1)*lda, 1 );
+ }
+
}
-#ifdef PRINTOUT
- IF (JJ==II)POUT(N7,J,K,X,X,Y,Y);
-#endif
-
-/*
-********************************************
-* Module 8: Procedure calls
-********************************************
-*/
- X = 1.0;
- Y = 1.0;
- Z = 1.0;
-
- for (I = 1; I <= N8; I++)
- {
- P3(X,Y,&Z);
- }
-#ifdef PRINTOUT
- IF (JJ==II)POUT(N8,J,K,X,Y,Z,Z);
-#endif
+ ipvt[n-1] = n;
-/*
-********************************************
-* Module 9: Array references
-********************************************
-*/
- J = 1;
- K = 2;
- L = 3;
- E1[1] = 1.0;
- E1[2] = 2.0;
- E1[3] = 3.0;
-
- for (I = 1; I <= N9; I++)
- {
- P0();
- }
-#ifdef PRINTOUT
- IF (JJ==II) POUT(N9,J,K,E1[1],E1[2],E1[3],E1[4]);
-#endif
-
-/*
-********************************************
-* Module 10: Integer arithmetic
-********************************************
-*/
- J = 2;
- K = 3;
-
- for (I = 1; I <= N10; I++)
- {
- J = J + K;
- K = J + K;
- J = K - J;
- K = K - J - J;
+ if ( a[n-1+(n-1)*lda] == 0.0 ) {
+ info = n;
}
-#ifdef PRINTOUT
- IF (JJ==II) POUT(N10,J,K,X1,X2,X3,X4);
-#endif
+ return info;
+}
+/******************************************************************************/
+void dgesl ( double a[], int lda, int n, int ipvt[], double b[], int job )
+
+/******************************************************************************/
/*
-********************************************
-* Module 11: Standard functions
-********************************************
-*/
- X = 0.75;
+ Purpose:
+
+ DGESL solves a real general linear system A * X = B.
+
+ Discussion:
+
+ DGESL can solve either of the systems A * X = B or A' * X = B.
+
+ The system matrix must have been factored by DGECO or DGEFA.
+
+ A division by zero will occur if the input factor contains a
+ zero on the diagonal. Technically this indicates singularity
+ but it is often caused by improper arguments or improper
+ setting of LDA. It will not occur if the subroutines are
+ called correctly and if DGECO has set 0.0 < RCOND
+ or DGEFA has set INFO == 0.
+
+ Modified:
+
+ 16 May 2005
+
+ Author:
+
+ C version by John Burkardt.
+
+ Reference:
+
+ Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
+ LINPACK User's Guide,
+ SIAM, (Society for Industrial and Applied Mathematics),
+ 3600 University City Science Center,
+ Philadelphia, PA, 19104-2688.
+ ISBN 0-89871-172-X
+
+ Parameters:
+
+ Input, double A[LDA*N], the output from DGECO or DGEFA.
+
+ Input, int LDA, the leading dimension of A.
+
+ Input, int N, the order of the matrix A.
+
+ Input, int IPVT[N], the pivot vector from DGECO or DGEFA.
+
+ Input/output, double B[N].
+ On input, the right hand side vector.
+ On output, the solution vector.
- for (I = 1; I <= N11; I++)
- {
- X = DSQRT(DEXP(DLOG(X)/T1));
- }
-#ifdef PRINTOUT
- IF (JJ==II) POUT(N11,J,K,X,X,X,X);
-#endif
+ Input, int JOB.
+ 0, solve A * X = B;
+ nonzero, solve A' * X = B.
+*/
+{
+ int k;
+ int l;
+ double t;
+ /*
+ Solve A * X = B.
+ */
+ if ( job == 0 ) {
+ for ( k = 1; k <= n-1; k++ ) {
+ l = ipvt[k-1];
+ t = b[l-1];
+
+ if ( l != k ) {
+ b[l-1] = b[k-1];
+ b[k-1] = t;
+ }
+
+ daxpy ( n-k, t, a+k+(k-1)*lda, 1, b+k, 1 );
+
+ }
-/*
-********************************************
-* THIS IS THE END OF THE MAJOR LOOP.
-********************************************
-*/
- if (++JJ <= II)
- goto IILOOP;
+ for ( k = n; 1 <= k; k-- ) {
+ b[k-1] = b[k-1] / a[k-1+(k-1)*lda];
+ t = -b[k-1];
+ daxpy ( k-1, t, a+0+(k-1)*lda, 1, b, 1 );
+ }
+ }
+ /*
+ Solve A' * X = B.
+ */
+ else {
+ for ( k = 1; k <= n; k++ ) {
+ t = ddot ( k-1, a+0+(k-1)*lda, 1, b, 1 );
+ b[k-1] = ( b[k-1] - t ) / a[k-1+(k-1)*lda];
+ }
+ for ( k = n-1; 1 <= k; k-- ) {
+ b[k-1] = b[k-1] + ddot ( n-k, a+k+(k-1)*lda, 1, b+k, 1 );
+ l = ipvt[k-1];
+
+ if ( l != k ) {
+ t = b[l-1];
+ b[l-1] = b[k-1];
+ b[k-1] = t;
+ }
+ }
+ }
+ return;
+}
+/******************************************************************************/
+
+void dscal ( int n, double sa, double x[], int incx )
+
+/******************************************************************************/
/*
-********************************************
-* Stop benchmark timing at this point.
-********************************************
+ Purpose:
+
+ DSCAL scales a vector by a constant.
+
+ Modified:
+
+ 30 March 2007
+
+ Author:
+
+ FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
+ C version by John Burkardt
+
+ Reference:
+
+ Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
+ LINPACK User's Guide,
+ SIAM, 1979.
+
+ Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
+ Basic Linear Algebra Subprograms for Fortran Usage,
+ Algorithm 539,
+ ACM Transactions on Mathematical Software,
+ Volume 5, Number 3, September 1979, pages 308-323.
+
+ Parameters:
+
+ Input, int N, the number of entries in the vector.
+
+ Input, double SA, the multiplier.
+
+ Input/output, double X[*], the vector to be scaled.
+
+ Input, int INCX, the increment between successive entries of X.
*/
- // finisec = time(0);
- finisec = timer.read_us();
- //timer.reset();
+{
+ int i;
+ int ix;
+ int m;
+
+ if ( n <= 0 ) {
+ } else if ( incx == 1 ) {
+ m = n % 5;
+
+ for ( i = 0; i < m; i++ ) {
+ x[i] = sa * x[i];
+ }
+ for ( i = m; i < n; i = i + 5 ) {
+ x[i] = sa * x[i];
+ x[i+1] = sa * x[i+1];
+ x[i+2] = sa * x[i+2];
+ x[i+3] = sa * x[i+3];
+ x[i+4] = sa * x[i+4];
+ }
+ } else {
+ if ( 0 <= incx ) {
+ ix = 0;
+ } else {
+ ix = ( - n + 1 ) * incx;
+ }
+
+ for ( i = 0; i < n; i++ ) {
+ x[ix] = sa * x[ix];
+ ix = ix + incx;
+ }
+ }
+ return;
+}
+/******************************************************************************/
+
+int idamax ( int n, double dx[], int incx )
+
+/******************************************************************************/
/*
-*--------------------------------------------------------------------
-* Performance in Whetstone KIP's per second is given by
-*
-* (100*LOOP*II)/TIME
-*
-* where TIME is in seconds.
-*--------------------------------------------------------------------
+ Purpose:
+
+ IDAMAX finds the index of the vector element of maximum absolute value.
+
+ Discussion:
+
+ WARNING: This index is a 1-based index, not a 0-based index!
+
+ Modified:
+
+ 30 March 2007
+
+ Author:
+
+ FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
+ C version by John Burkardt
+
+ Reference:
+
+ Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
+ LINPACK User's Guide,
+ SIAM, 1979.
+
+ Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
+ Basic Linear Algebra Subprograms for Fortran Usage,
+ Algorithm 539,
+ ACM Transactions on Mathematical Software,
+ Volume 5, Number 3, September 1979, pages 308-323.
+
+ Parameters:
+
+ Input, int N, the number of entries in the vector.
+
+ Input, double X[*], the vector to be examined.
+
+ Input, int INCX, the increment between successive entries of SX.
+
+ Output, int IDAMAX, the index of the element of maximum
+ absolute value.
*/
- pc.printf(" kraj \n");
- double vreme;
- vreme = (finisec - startsec) / 1000000;
-
- if (vreme <= 0)
- {
- pc.printf("Insufficient duration- Increase the LOOP count \n");
- finisec = 0;
- startsec = 0;
- return 1;
- }
+{
+ double dmax;
+ int i;
+ int ix;
+ int value;
+
+ value = 0;
- pc.printf("Loops: %ld , \t Iterations: %d, \t Duration: %.3f sec. \n",
- LOOP, II, vreme);
+ if ( n < 1 || incx <= 0 ) {
+ return value;
+ }
+
+ value = 1;
+
+ if ( n == 1 ) {
+ return value;
+ }
- KIPS = (100.0 * LOOP * II) / vreme ;
-
- // if (KIPS >= 1000.0)
- // pc.printf("C Converted Double Precision Whetstones: %.3f MIPS \n\n", KIPS / 1000);
- // else
- // pc.printf("C Converted Double Precision Whetstones: %.3f KIPS \n\n", KIPS);
-
- pc.printf("C Converted Double Precision Whetstones: %.3f MIPS \n\n", KIPS / 1000);
+ if ( incx == 1 ) {
+ dmax = r8_abs ( dx[0] );
+
+ for ( i = 1; i < n; i++ ) {
+ if ( dmax < r8_abs ( dx[i] ) ) {
+ value = i + 1;
+ dmax = r8_abs ( dx[i] );
+ }
+ }
+ } else {
+ ix = 0;
+ dmax = r8_abs ( dx[0] );
+ ix = ix + incx;
+
+ for ( i = 1; i < n; i++ ) {
+ if ( dmax < r8_abs ( dx[ix] ) ) {
+ value = i + 1;
+ dmax = r8_abs ( dx[ix] );
+ }
+ ix = ix + incx;
+ }
+ }
- if (continuous)
- goto LCONT;
+ return value;
+}
+/******************************************************************************/
+
+double r8_abs ( double x )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ R8_ABS returns the absolute value of a R8.
- finisec = 0;
- startsec = 0;
- return 1;
-}
+ Modified:
+
+ 02 April 2005
+
+ Author:
-void PA(double E[])
+ John Burkardt
+
+ Parameters:
+
+ Input, double X, the quantity whose absolute value is desired.
+
+ Output, double R8_ABS, the absolute value of X.
+*/
{
- J = 0;
+ double value;
+
+ if ( 0.0 <= x ) {
+ value = x;
+ } else {
+ value = -x;
+ }
+ return value;
+}
+/******************************************************************************/
+
+double r8_epsilon ( void )
-L10:
- E[1] = ( E[1] + E[2] + E[3] - E[4]) * T;
- E[2] = ( E[1] + E[2] - E[3] + E[4]) * T;
- E[3] = ( E[1] - E[2] + E[3] + E[4]) * T;
- E[4] = (-E[1] + E[2] + E[3] + E[4]) / T2;
- J += 1;
+/******************************************************************************/
+/*
+ Purpose:
+
+ R8_EPSILON returns the R8 round off unit.
+
+ Discussion:
+
+ R8_EPSILON is a number R which is a power of 2 with the property that,
+ to the precision of the computer's arithmetic,
+ 1 < 1 + R
+ but
+ 1 = ( 1 + R / 2 )
- if (J < 6)
- goto L10;
-}
+ Licensing:
+
+ This code is distributed under the GNU LGPL license.
+
+ Modified:
+
+ 08 May 2006
+
+ Author:
+
+ John Burkardt
+
+ Parameters:
-void P0(void)
+ Output, double R8_EPSILON, the double precision round-off unit.
+*/
{
- E1[J] = E1[K];
- E1[K] = E1[L];
- E1[L] = E1[J];
+ double r;
+
+ r = 1.0;
+
+ while ( 1.0 < ( double ) ( 1.0 + r ) ) {
+ r = r / 2.0;
+ }
+ r = 2.0 * r;
+
+ return r;
}
+/******************************************************************************/
-void P3(double X, double Y, double *Z)
-{
- double X1, Y1;
+double r8_max ( double x, double y )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ R8_MAX returns the maximum of two R8's.
+
+ Modified:
+
+ 18 August 2004
+
+ Author:
- X1 = X;
- Y1 = Y;
- X1 = T * (X1 + Y1);
- Y1 = T * (X1 + Y1);
- *Z = (X1 + Y1) / T2;
-}
+ John Burkardt
+
+ Parameters:
+
+ Input, double X, Y, the quantities to compare.
-#ifdef PRINTOUT
-void POUT(long N, long J, long K, double X1, double X2, double X3, double X4)
+ Output, double R8_MAX, the maximum of X and Y.
+*/
{
- pc.printf("%7ld %7ld %7ld %12.4e %12.4e %12.4e %12.4e\n",
- N, J, K, X1, X2, X3, X4);
+ double value;
+
+ if ( y < x ) {
+ value = x;
+ } else {
+ value = y;
+ }
+ return value;
}
-#endif
+/******************************************************************************/
+
+double r8_random ( int iseed[4] )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ R8_RANDOM returns a uniformly distributed random number between 0 and 1.
+
+ Discussion:
-//*********************************
-//** MAIN block **
-//*********************************
-int main()
+ This routine uses a multiplicative congruential method with modulus
+ 2**48 and multiplier 33952834046453 (see G.S.Fishman,
+ 'Multiplicative congruential random number generators with modulus
+ 2**b: an exhaustive analysis for b = 32 and a partial analysis for
+ b = 48', Math. Comp. 189, pp 331-344, 1990).
+
+ 48-bit integers are stored in 4 integer array elements with 12 bits
+ per element. Hence the routine is portable across machines with
+ integers of 32 bits or more.
+
+ Parameters:
+
+ Input/output, integer ISEED(4).
+ On entry, the seed of the random number generator; the array
+ elements must be between 0 and 4095, and ISEED(4) must be odd.
+ On exit, the seed is updated.
+
+ Output, double R8_RANDOM, the next pseudorandom number.
+*/
{
- int rez=0;
- printf("\n My Benchamrk example for Whetstones \n");
-
- while(1)
- {
- myled=1-rez;
-
- rez = Whetstone(); //Call of Whetstone banch metod
-
- myled=1-rez;
- wait_us(0.3);
+ int ipw2 = 4096;
+ int it1;
+ int it2;
+ int it3;
+ int it4;
+ int m1 = 494;
+ int m2 = 322;
+ int m3 = 2508;
+ int m4 = 2549;
+ double r = 1.0 / 4096.0;
+ double value;
+ /*
+ Multiply the seed by the multiplier modulo 2**48.
+ */
+ it4 = iseed[3] * m4;
+ it3 = it4 / ipw2;
+ it4 = it4 - ipw2 * it3;
+ it3 = it3 + iseed[2] * m4 + iseed[3] * m3;
+ it2 = it3 / ipw2;
+ it3 = it3 - ipw2 * it2;
+ it2 = it2 + iseed[1] * m4 + iseed[2] * m3 + iseed[3] * m2;
+ it1 = it2 / ipw2;
+ it2 = it2 - ipw2 * it1;
+ it1 = it1 + iseed[0] * m4 + iseed[1] * m3 + iseed[2] * m2 + iseed[3] * m1;
+ it1 = ( it1 % ipw2 );
+ /*
+ Return updated seed
+ */
+ iseed[0] = it1;
+ iseed[1] = it2;
+ iseed[2] = it3;
+ iseed[3] = it4;
+ /*
+ Convert 48-bit integer to a real number in the interval (0,1)
+ */
+ value =
+ r * ( ( double ) ( it1 )
+ + r * ( ( double ) ( it2 )
+ + r * ( ( double ) ( it3 )
+ + r * ( ( double ) ( it4 ) ) ) ) );
+
+ return value;
+}
+/******************************************************************************/
+
+double *r8mat_gen ( int lda, int n )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ R8MAT_GEN generates a random R8MAT.
+
+ Modified:
+
+ 06 June 2005
+
+ Parameters:
+
+ Input, integer LDA, the leading dimension of the matrix.
+
+ Input, integer N, the order of the matrix.
+
+ Output, double R8MAT_GEN[LDA*N], the N by N matrix.
+*/
+{
+ double *a;
+ int i;
+ int init[4] = { 1, 2, 3, 1325 };
+ int j;
+
+ a = ( double * ) malloc ( lda * n * sizeof ( double ) );
+
+ for ( j = 1; j <= n; j++ ) {
+ for ( i = 1; i <= n; i++ ) {
+ a[i-1+(j-1)*lda] = r8_random ( init ) - 0.5;
+ }
}
+
+ return a;
}
\ No newline at end of file