Jovan Ivković / Mbed 2 deprecated 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.

Dependencies:   mbed

main.cpp@2:273153e44338, 2017-01-03 (annotated)

Committer:
JovanEps
Date:
Tue Jan 03 02:31:38 2017 +0000
Revision:
2:273153e44338
Parent:
1:be78b18b8347
Child:
3:da1132c65314
1v

Who changed what in which revision?

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