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@10:7473e1529a9e, 2017-04-09 (annotated)

Committer:
JovanEps
Date:
Sun Apr 09 20:31:20 2017 +0000
Revision:
10:7473e1529a9e
Parent:
9:54628dc6805e 
B3

Who changed what in which revision?

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