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@6:5e0f3eedaf66, 2017-01-03 (annotated)

Committer:
JovanEps
Date:
Tue Jan 03 19:59:24 2017 +0000
Revision:
6:5e0f3eedaf66
Parent:
5:2b929fbd5c69
Child:
7:931219974070
RC2

Who changed what in which revision?

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