← Back to branch summary

~ubuntu-branches/ubuntu/vivid/emscripten/vivid

« back to all changes in this revision

Viewing changes to tests/linpack.c

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2013-06-11 15:45:24 UTC
  • mfrom: (1.2.1) (2.1.1 experimental)
  • Revision ID: package-import@ubuntu.com-20130611154524-rppb3w6tixlegv4n
Tags: 1.4.7~20130611~a1eb425-1
* New snapshot release
* Upload to unstable

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/linpack.c
 
1
/*      gcc linpack.c cpuidc64.o cpuida64.o -m64 -lrt -lc -lm -o linpack
 
2
 *
 
3
 *          Linpack 100x100 Benchmark In C/C++ For PCs
 
4
 *     
 
5
 *  Different compilers can produce different floating point numeric
 
6
 *  results, probably due to compiling instructions in a different
 
7
 *  sequence. As the program checks these, they may need to be changed.
 
8
 *  The log file indicates non-standard results and these values can
 
9
 *  be copied and pasted into this program. See // Values near the
 
10
 *  end of main().
 
11
 *
 
12
 *  Different compilers do not optimise the code in the same way.
 
13
 *  This can lead to wide variations in benchmark speeds. See results
 
14
 *  with MS6 compiler ID and compare with those from same CPUs from
 
15
 *  the Watcom compiler generated code.
 
16
 *
 
17
 ***************************************************************************
 
18
*/
 
19
 
 
20
#define _CRT_SECURE_NO_WARNINGS 1
 
21
#ifdef WIN32
 
22
#include <Windows.h>
 
23
#else
 
24
#include <sys/time.h>
 
25
#endif
 
26
 
 
27
#define UNROLL
 
28
#define DP
 
29
 
 
30
#ifdef SP
 
31
#define REAL float
 
32
#define ZERO 0.0
 
33
#define ONE 1.0
 
34
#define PREC "Single"
 
35
#endif
 
36
 
 
37
#ifdef DP
 
38
#define REAL double
 
39
#define ZERO 0.0e0
 
40
#define ONE 1.0e0
 
41
#define PREC "Double"
 
42
#endif
 
43
 
 
44
#ifdef ROLL
 
45
#define ROLLING "Rolled"
 
46
#endif
 
47
#ifdef UNROLL
 
48
#define ROLLING "Unrolled"
 
49
#endif
 
50
 
 
51
 // VERSION
 
52
               
 
53
 #ifdef CNNT
 
54
    #define options   "Non-optimised"
 
55
    #define opt "0"
 
56
 #else
 
57
//    #define options   "Optimised"
 
58
    #define options   "Opt 3 64 Bit"
 
59
    #define opt "1"
 
60
 #endif
 
61
 
 
62
#define NTIMES 10
 
63
 
 
64
#include <stdio.h>
 
65
#include <math.h>
 
66
#include <stdlib.h>
 
67
#include <time.h> 
 
68
 
 
69
 
 
70
/* this is truly rank, but it's minimally invasive, and lifted in part from the STREAM scores */
 
71
 
 
72
static double secs;
 
73
 
 
74
#ifndef WIN32
 
75
 
 
76
double mysecond()
 
77
{
 
78
        struct timeval tp;
 
79
        struct timezone tzp;
 
80
        int i;
 
81
 
 
82
        i = gettimeofday(&tp,&tzp);
 
83
        return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 );
 
84
}
 
85
#else
 
86
 
 
87
double mysecond()
 
88
{
 
89
        static LARGE_INTEGER freq = {0};
 
90
        LARGE_INTEGER count = {0};
 
91
        if(freq.QuadPart == 0LL) {
 
92
                QueryPerformanceFrequency(&freq);
 
93
        }
 
94
        QueryPerformanceCounter(&count);
 
95
        return (double)count.QuadPart / (double)freq.QuadPart;
 
96
}
 
97
 
 
98
#endif
 
99
 
 
100
void start_time()
 
101
{
 
102
        secs = mysecond();
 
103
}
 
104
 
 
105
void end_time()
 
106
{
 
107
        secs = mysecond() - secs;
 
108
}
 
109
 
 
110
void print_time (int row);
 
111
void matgen (REAL a[], int lda, int n, REAL b[], REAL *norma);
 
112
void dgefa (REAL a[], int lda, int n, int ipvt[], int *info);
 
113
void dgesl (REAL a[],int lda,int n,int ipvt[],REAL b[],int job);
 
114
void dmxpy (int n1, REAL y[], int n2, int ldm, REAL x[], REAL m[]);
 
115
void daxpy (int n, REAL da, REAL dx[], int incx, REAL dy[], int incy);
 
116
REAL epslon (REAL x);
 
117
int idamax (int n, REAL dx[], int incx);
 
118
void dscal (int n, REAL da, REAL dx[], int incx);
 
119
REAL ddot (int n, REAL dx[], int incx, REAL dy[], int incy);
 
120
 
 
121
static REAL atime[9][15];
 
122
double runSecs = 1;
 
123
 
 
124
 
 
125
int main (int argc, char *argv[])
 
126
{
 
127
        static REAL aa[200*200],a[200*201],b[200],x[200];       
 
128
        REAL cray,ops,total,norma,normx;
 
129
        REAL resid,residn,eps,tm2,epsn,x1,x2;
 
130
        REAL mflops;
 
131
        static int ipvt[200],n,i,j,ntimes,info,lda,ldaa;
 
132
        int endit, pass, loop;
 
133
        REAL overhead1, overhead2, time2;
 
134
        REAL max1, max2;
 
135
        char was[5][20];
 
136
        char expect[5][20];
 
137
        char title[5][20];
 
138
        int errors;
 
139
        
 
140
 
 
141
        printf("\n");
 
142
         
 
143
        printf("##########################################\n"); 
 
144
 
 
145
 
 
146
    
 
147
        lda = 201;
 
148
        ldaa = 200;
 
149
        cray = .056; 
 
150
        n = 100;
 
151
 
 
152
        fprintf(stdout, "%s ", ROLLING);
 
153
        fprintf(stdout, "%s ", PREC);
 
154
        fprintf(stdout,"Precision Linpack Benchmark - PC Version in 'C/C++'\n\n");
 
155
 
 
156
        fprintf(stdout,"Optimisation %s\n\n",options);
 
157
 
 
158
        ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);
 
159
 
 
160
        matgen(a,lda,n,b,&norma);
 
161
        start_time();
 
162
        dgefa(a,lda,n,ipvt,&info);
 
163
        end_time();
 
164
        atime[0][0] = secs;
 
165
        start_time();
 
166
        dgesl(a,lda,n,ipvt,b,0);
 
167
        end_time();
 
168
        atime[1][0] = secs;
 
169
        total = atime[0][0] + atime[1][0];
 
170
 
 
171
/*     compute a residual to verify results.  */ 
 
172
 
 
173
        for (i = 0; i < n; i++) {
 
174
                x[i] = b[i];
 
175
        }
 
176
        matgen(a,lda,n,b,&norma);
 
177
        for (i = 0; i < n; i++) {
 
178
                b[i] = -b[i];
 
179
        }
 
180
        dmxpy(n,b,n,lda,x,a);
 
181
        resid = 0.0;
 
182
        normx = 0.0;
 
183
        for (i = 0; i < n; i++) {
 
184
                resid = (resid > fabs((double)b[i])) 
 
185
                        ? resid : fabs((double)b[i]);
 
186
                normx = (normx > fabs((double)x[i])) 
 
187
                        ? normx : fabs((double)x[i]);
 
188
        }
 
189
        eps = epslon(ONE);
 
190
        residn = resid/( n*norma*normx*eps );
 
191
        epsn = eps;
 
192
        x1 = x[0] - 1;
 
193
        x2 = x[n-1] - 1;
 
194
        
 
195
        printf("norm resid      resid           machep");
 
196
        printf("         x[0]-1          x[n-1]-1\n");
 
197
        printf("%6.1f %17.8e%17.8e%17.8e%17.8e\n\n",
 
198
               (double)residn, (double)resid, (double)epsn, 
 
199
               (double)x1, (double)x2);
 
200
 
 
201
        printf("Times are reported for matrices of order        %5d\n",n);
 
202
        printf("1 pass times for array with leading dimension of%5d\n\n",lda);
 
203
        printf("      dgefa      dgesl      total     Mflops       unit");
 
204
        printf("      ratio\n");
 
205
 
 
206
        atime[2][0] = total;
 
207
        if (total > 0.0)
 
208
        {
 
209
            atime[3][0] = ops/(1.0e6*total);
 
210
            atime[4][0] = 2.0/atime[3][0];
 
211
        }
 
212
        else
 
213
        {
 
214
            atime[3][0] = 0.0;
 
215
            atime[4][0] = 0.0;
 
216
        }
 
217
        atime[5][0] = total/cray;
 
218
       
 
219
        print_time(0);
 
220
 
 
221
/************************************************************************
 
222
 *       Calculate overhead of executing matgen procedure              *
 
223
 ************************************************************************/
 
224
       
 
225
        printf("\nCalculating matgen overhead\n");
 
226
        pass = -20;
 
227
        loop = NTIMES;
 
228
        do
 
229
        {
 
230
            start_time();
 
231
            pass = pass + 1;        
 
232
            for ( i = 0 ; i < loop ; i++)
 
233
            {
 
234
                 matgen(a,lda,n,b,&norma);
 
235
            }
 
236
            end_time();
 
237
            overhead1 = secs;
 
238
            printf("%10d times %6.2f seconds\n", loop, overhead1);
 
239
            if (overhead1 > runSecs)
 
240
            {
 
241
                pass = 0;
 
242
            }
 
243
            if (pass < 0)
 
244
            {
 
245
                if (overhead1 < 0.1)
 
246
                {
 
247
                    loop = loop * 10;
 
248
                }
 
249
                else
 
250
                {
 
251
                    loop = loop * 2;
 
252
                }
 
253
            }
 
254
        }
 
255
        while (pass < 0);
 
256
        
 
257
        overhead1 = overhead1 / (double)loop;
 
258
 
 
259
        printf("Overhead for 1 matgen %12.5f seconds\n\n", overhead1);
 
260
 
 
261
/************************************************************************
 
262
 *           Calculate matgen/dgefa passes for runSecs seconds                *
 
263
 ************************************************************************/
 
264
       
 
265
        printf("Calculating matgen/dgefa passes for %d seconds\n", (int)runSecs);
 
266
        pass = -20;
 
267
        ntimes = NTIMES;
 
268
        do
 
269
        {
 
270
            start_time();
 
271
            pass = pass + 1;        
 
272
            for ( i = 0 ; i < ntimes ; i++)
 
273
            {
 
274
                matgen(a,lda,n,b,&norma);
 
275
                dgefa(a,lda,n,ipvt,&info );
 
276
            }
 
277
            end_time();
 
278
            time2 = secs;
 
279
            printf("%10d times %6.2f seconds\n", ntimes, time2);
 
280
            if (time2 > runSecs)
 
281
            {
 
282
                pass = 0;
 
283
            }
 
284
            if (pass < 0)
 
285
            {
 
286
                if (time2 < 0.1)
 
287
                {
 
288
                    ntimes = ntimes * 10;
 
289
                }
 
290
                else
 
291
                {
 
292
                    ntimes = ntimes * 2;
 
293
                }
 
294
            }
 
295
        }
 
296
        while (pass < 0);
 
297
        
 
298
        ntimes =  (int)(runSecs * (double)ntimes / time2);
 
299
        if (ntimes == 0) ntimes = 1;
 
300
 
 
301
        printf("Passes used %10d \n\n", ntimes);
 
302
        printf("Times for array with leading dimension of%4d\n\n",lda);
 
303
        printf("      dgefa      dgesl      total     Mflops       unit");
 
304
        printf("      ratio\n");        
 
305
 
 
306
/************************************************************************
 
307
 *                              Execute 5 passes                        *
 
308
 ************************************************************************/
 
309
      
 
310
        tm2 = ntimes * overhead1;
 
311
        atime[3][6] = 0;
 
312
 
 
313
        for (j=1 ; j<6 ; j++)
 
314
        {
 
315
            start_time();
 
316
            for (i = 0; i < ntimes; i++)
 
317
            {
 
318
                matgen(a,lda,n,b,&norma);
 
319
                dgefa(a,lda,n,ipvt,&info );
 
320
            }
 
321
            end_time();
 
322
            atime[0][j] = (secs - tm2)/ntimes;
 
323
 
 
324
            start_time();              
 
325
            for (i = 0; i < ntimes; i++)
 
326
            {
 
327
                dgesl(a,lda,n,ipvt,b,0);
 
328
            }
 
329
            end_time();
 
330
 
 
331
            atime[1][j] = secs/ntimes;
 
332
            total       = atime[0][j] + atime[1][j];
 
333
            atime[2][j] = total;
 
334
            atime[3][j] = ops/(1.0e6*total);
 
335
            atime[4][j] = 2.0/atime[3][j];
 
336
            atime[5][j] = total/cray;
 
337
            atime[3][6] = atime[3][6] + atime[3][j];
 
338
            
 
339
            print_time(j);
 
340
        }
 
341
        atime[3][6] = atime[3][6] / 5.0;
 
342
        printf("Average                          %11.2f\n",
 
343
                                               (double)atime[3][6]);        
 
344
        
 
345
        printf("\nCalculating matgen2 overhead\n");
 
346
 
 
347
/************************************************************************
 
348
 *             Calculate overhead of executing matgen procedure         *
 
349
 ************************************************************************/
 
350
 
 
351
        start_time();        
 
352
        for ( i = 0 ; i < loop ; i++)
 
353
        {
 
354
            matgen(aa,ldaa,n,b,&norma);    
 
355
        }
 
356
        end_time();
 
357
        overhead2 = secs;
 
358
        overhead2 = overhead2 / (double)loop;
 
359
        
 
360
        printf("Overhead for 1 matgen %12.5f seconds\n\n", overhead2);
 
361
        printf("Times for array with leading dimension of%4d\n\n",ldaa);
 
362
        printf("      dgefa      dgesl      total     Mflops       unit");
 
363
        printf("      ratio\n");
 
364
 
 
365
/************************************************************************
 
366
 *                              Execute 5 passes                        *
 
367
 ************************************************************************/
 
368
              
 
369
        tm2 = ntimes * overhead2;
 
370
        atime[3][12] = 0;
 
371
 
 
372
        for (j=7 ; j<12 ; j++)
 
373
        {
 
374
            start_time();
 
375
            for (i = 0; i < ntimes; i++)
 
376
            {
 
377
                matgen(aa,ldaa,n,b,&norma);
 
378
                dgefa(aa,ldaa,n,ipvt,&info  );
 
379
            }
 
380
            end_time();
 
381
            atime[0][j] = (secs - tm2)/ntimes;
 
382
            
 
383
            start_time();      
 
384
            for (i = 0; i < ntimes; i++)
 
385
            {
 
386
                dgesl(aa,ldaa,n,ipvt,b,0);
 
387
            }
 
388
            end_time();
 
389
            atime[1][j] = secs/ntimes;
 
390
            total       = atime[0][j] + atime[1][j];
 
391
            atime[2][j] = total;
 
392
            atime[3][j] = ops/(1.0e6*total);
 
393
            atime[4][j] = 2.0/atime[3][j];
 
394
            atime[5][j] = total/cray;
 
395
            atime[3][12] = atime[3][12] + atime[3][j];
 
396
 
 
397
            print_time(j);
 
398
        }
 
399
        atime[3][12] = atime[3][12] / 5.0; 
 
400
        printf("Average                          %11.2f\n",
 
401
                                              (double)atime[3][12]);  
 
402
 
 
403
/************************************************************************
 
404
 *           Use minimum average as overall Mflops rating               *
 
405
 ************************************************************************/
 
406
      
 
407
        mflops = atime[3][6];
 
408
        if (atime[3][12] < mflops) mflops = atime[3][12];
 
409
       
 
410
        printf("\n");
 
411
        printf( "%s ", ROLLING);
 
412
        printf( "%s ", PREC);
 
413
        printf(" Precision %11.2f Mflops \n\n",mflops);
 
414
 
 
415
 
 
416
    max1 = 0;
 
417
    for (i=1 ; i<6 ; i++)
 
418
    {
 
419
        if (atime[3][i] > max1) max1 = atime[3][i];                 
 
420
    }
 
421
 
 
422
    max2 = 0;
 
423
    for (i=7 ; i<12 ; i++)
 
424
    {                 
 
425
        if (atime[3][i] > max2) max2 = atime[3][i];                 
 
426
    }
 
427
    if (max1 < max2) max2 = max1;
 
428
   
 
429
    sprintf(was[0], "%16.1f",(double)residn);
 
430
    sprintf(was[1], "%16.8e",(double)resid);
 
431
    sprintf(was[2], "%16.8e",(double)epsn);
 
432
    sprintf(was[3], "%16.8e",(double)x1);
 
433
    sprintf(was[4], "%16.8e",(double)x2);
 
434
 
 
435
/*
 
436
    //  Values for Watcom
 
437
 
 
438
    sprintf(expect[0], "             0.4");
 
439
    sprintf(expect[1], " 7.41628980e-014");
 
440
    sprintf(expect[2], " 1.00000000e-015");
 
441
    sprintf(expect[3], "-1.49880108e-014");
 
442
    sprintf(expect[4], "-1.89848137e-014");
 
443
    // Values for Visual C++
 
444
 
 
445
    sprintf(expect[0], "             1.7");
 
446
    sprintf(expect[1], " 7.41628980e-014");
 
447
    sprintf(expect[2], " 2.22044605e-016");
 
448
    sprintf(expect[3], "-1.49880108e-014");
 
449
    sprintf(expect[4], "-1.89848137e-014");
 
450
 
 
451
    // Values for Ubuntu GCC 32 Bit
 
452
 
 
453
    sprintf(expect[0], "             1.9");
 
454
    sprintf(expect[1], "  8.39915160e-14");
 
455
    sprintf(expect[2], "  2.22044605e-16");
 
456
    sprintf(expect[3], " -6.22835117e-14");
 
457
    sprintf(expect[4], " -4.16333634e-14");
 
458
*/
 
459
 
 
460
     // Values for Ubuntu GCC 32 Bit
 
461
 
 
462
    sprintf(expect[0], "             1.7");
 
463
    sprintf(expect[1], "  7.41628980e-14");
 
464
    sprintf(expect[2], "  2.22044605e-16");
 
465
    sprintf(expect[3], " -1.49880108e-14");
 
466
    sprintf(expect[4], " -1.89848137e-14");
 
467
 
 
468
    sprintf(title[0], "norm. resid");
 
469
    sprintf(title[1], "resid      ");
 
470
    sprintf(title[2], "machep     ");
 
471
    sprintf(title[3], "x[0]-1     ");
 
472
    sprintf(title[4], "x[n-1]-1   ");
 
473
 
 
474
    if (strtol(opt, NULL, 10) == 0)
 
475
    {
 
476
        sprintf(expect[2], " 8.88178420e-016");
 
477
    }
 
478
    errors = 0;
 
479
 
 
480
    printf ("\n");
 
481
}
 
482
     
 
483
/*----------------------*/ 
 
484
void print_time (int row)
 
485
 
 
486
{
 
487
printf("%11.5f%11.5f%11.5f%11.2f%11.4f%11.4f\n",   (double)atime[0][row],
 
488
       (double)atime[1][row], (double)atime[2][row], (double)atime[3][row], 
 
489
       (double)atime[4][row], (double)atime[5][row]);
 
490
       return;
 
491
}
 
492
      
 
493
/*----------------------*/ 
 
494
 
 
495
void matgen (REAL a[], int lda, int n, REAL b[], REAL *norma)
 
496
 
 
497
 
 
498
/* We would like to declare a[][lda], but c does not allow it.  In this
 
499
function, references to a[i][j] are written a[lda*i+j].  */
 
500
 
 
501
{
 
502
        int init, i, j;
 
503
 
 
504
        init = 1325;
 
505
        *norma = 0.0;
 
506
        for (j = 0; j < n; j++) {
 
507
                for (i = 0; i < n; i++) {
 
508
                        init = 3125*init % 65536;
 
509
                        a[lda*j+i] = (init - 32768.0)/16384.0;                        
 
510
                        *norma = (a[lda*j+i] > *norma) ? a[lda*j+i] : *norma;
 
511
                        
 
512
                        /* alternative for some compilers
 
513
                        if (fabs(a[lda*j+i]) > *norma) *norma = fabs(a[lda*j+i]);
 
514
                        */
 
515
                }
 
516
        }
 
517
        for (i = 0; i < n; i++) {
 
518
          b[i] = 0.0;
 
519
        }
 
520
        for (j = 0; j < n; j++) {
 
521
                for (i = 0; i < n; i++) {
 
522
                        b[i] = b[i] + a[lda*j+i];
 
523
                }
 
524
        }
 
525
        return;
 
526
}
 
527
 
 
528
/*----------------------*/ 
 
529
void dgefa(REAL a[], int lda, int n, int ipvt[], int *info)
 
530
 
 
531
 
 
532
/* We would like to declare a[][lda], but c does not allow it.  In this
 
533
function, references to a[i][j] are written a[lda*i+j].  */
 
534
/*
 
535
     dgefa factors a double precision matrix by gaussian elimination.
 
536
 
 
537
     dgefa is usually called by dgeco, but it can be called
 
538
     directly with a saving in time if  rcond  is not needed.
 
539
     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
 
540
 
 
541
     on entry
 
542
 
 
543
        a       REAL precision[n][lda]
 
544
                the matrix to be factored.
 
545
 
 
546
        lda     integer
 
547
                the leading dimension of the array  a .
 
548
 
 
549
        n       integer
 
550
                the order of the matrix  a .
 
551
 
 
552
     on return
 
553
 
 
554
        a       an upper triangular matrix and the multipliers
 
555
                which were used to obtain it.
 
556
                the factorization can be written  a = l*u  where
 
557
                l  is a product of permutation and unit lower
 
558
                triangular matrices and  u  is upper triangular.
 
559
 
 
560
        ipvt    integer[n]
 
561
                an integer vector of pivot indices.
 
562
 
 
563
        info    integer
 
564
                = 0  normal value.
 
565
                = k  if  u[k][k] .eq. 0.0 .  this is not an error
 
566
                     condition for this subroutine, but it does
 
567
                     indicate that dgesl or dgedi will divide by zero
 
568
                     if called.  use  rcond  in dgeco for a reliable
 
569
                     indication of singularity.
 
570
 
 
571
     linpack. this version dated 08/14/78 .
 
572
     cleve moler, university of new mexico, argonne national lab.
 
573
 
 
574
     functions
 
575
 
 
576
     blas daxpy,dscal,idamax
 
577
*/
 
578
 
 
579
{
 
580
/*     internal variables       */
 
581
 
 
582
REAL t;
 
583
int j,k,kp1,l,nm1;
 
584
 
 
585
 
 
586
/*     gaussian elimination with partial pivoting       */
 
587
 
 
588
        *info = 0;
 
589
        nm1 = n - 1;
 
590
        if (nm1 >=  0) {
 
591
                for (k = 0; k < nm1; k++) {
 
592
                        kp1 = k + 1;
 
593
 
 
594
                        /* find l = pivot index */
 
595
 
 
596
                        l = idamax(n-k,&a[lda*k+k],1) + k;
 
597
                        ipvt[k] = l;
 
598
 
 
599
                        /* zero pivot implies this column already 
 
600
                           triangularized */
 
601
 
 
602
                        if (a[lda*k+l] != ZERO) {
 
603
 
 
604
                                /* interchange if necessary */
 
605
 
 
606
                                if (l != k) {
 
607
                                        t = a[lda*k+l];
 
608
                                        a[lda*k+l] = a[lda*k+k];
 
609
                                        a[lda*k+k] = t; 
 
610
                                }
 
611
 
 
612
                                /* compute multipliers */
 
613
 
 
614
                                t = -ONE/a[lda*k+k];
 
615
                                dscal(n-(k+1),t,&a[lda*k+k+1],1);
 
616
 
 
617
                                /* row elimination with column indexing */
 
618
 
 
619
                                for (j = kp1; j < n; j++) {
 
620
                                        t = a[lda*j+l];
 
621
                                        if (l != k) {
 
622
                                                a[lda*j+l] = a[lda*j+k];
 
623
                                                a[lda*j+k] = t;
 
624
                                        }
 
625
                                        daxpy(n-(k+1),t,&a[lda*k+k+1],1,
 
626
                                              &a[lda*j+k+1],1);
 
627
                                } 
 
628
                        }
 
629
                        else { 
 
630
                                *info = k;
 
631
                        }
 
632
                } 
 
633
        }
 
634
        ipvt[n-1] = n-1;
 
635
        if (a[lda*(n-1)+(n-1)] == ZERO) *info = n-1;
 
636
        return;
 
637
}
 
638
 
 
639
/*----------------------*/ 
 
640
 
 
641
void dgesl(REAL a[],int lda,int n,int ipvt[],REAL b[],int job )
 
642
 
 
643
 
 
644
/* We would like to declare a[][lda], but c does not allow it.  In this
 
645
function, references to a[i][j] are written a[lda*i+j].  */
 
646
 
 
647
/*
 
648
     dgesl solves the double precision system
 
649
     a * x = b  or  trans(a) * x = b
 
650
     using the factors computed by dgeco or dgefa.
 
651
 
 
652
     on entry
 
653
 
 
654
        a       double precision[n][lda]
 
655
                the output from dgeco or dgefa.
 
656
 
 
657
        lda     integer
 
658
                the leading dimension of the array  a .
 
659
 
 
660
        n       integer
 
661
                the order of the matrix  a .
 
662
 
 
663
        ipvt    integer[n]
 
664
                the pivot vector from dgeco or dgefa.
 
665
 
 
666
        b       double precision[n]
 
667
                the right hand side vector.
 
668
 
 
669
        job     integer
 
670
                = 0         to solve  a*x = b ,
 
671
                = nonzero   to solve  trans(a)*x = b  where
 
672
                            trans(a)  is the transpose.
 
673
 
 
674
    on return
 
675
 
 
676
        b       the solution vector  x .
 
677
 
 
678
     error condition
 
679
 
 
680
        a division by zero will occur if the input factor contains a
 
681
        zero on the diagonal.  technically this indicates singularity
 
682
        but it is often caused by improper arguments or improper
 
683
        setting of lda .  it will not occur if the subroutines are
 
684
        called correctly and if dgeco has set rcond .gt. 0.0
 
685
        or dgefa has set info .eq. 0 .
 
686
 
 
687
     to compute  inverse(a) * c  where  c  is a matrix
 
688
     with  p  columns
 
689
           dgeco(a,lda,n,ipvt,rcond,z)
 
690
           if (!rcond is too small){
 
691
                for (j=0,j<p,j++)
 
692
                        dgesl(a,lda,n,ipvt,c[j][0],0);
 
693
           }
 
694
 
 
695
     linpack. this version dated 08/14/78 .
 
696
     cleve moler, university of new mexico, argonne national lab.
 
697
 
 
698
     functions
 
699
 
 
700
     blas daxpy,ddot
 
701
*/
 
702
{
 
703
/*     internal variables       */
 
704
 
 
705
        REAL t;
 
706
        int k,kb,l,nm1;
 
707
 
 
708
        nm1 = n - 1;
 
709
        if (job == 0) {
 
710
 
 
711
                /* job = 0 , solve  a * x = b
 
712
                   first solve  l*y = b         */
 
713
 
 
714
                if (nm1 >= 1) {
 
715
                        for (k = 0; k < nm1; k++) {
 
716
                                l = ipvt[k];
 
717
                                t = b[l];
 
718
                                if (l != k){ 
 
719
                                        b[l] = b[k];
 
720
                                        b[k] = t;
 
721
                                }       
 
722
                                daxpy(n-(k+1),t,&a[lda*k+k+1],1,&b[k+1],1 );
 
723
                        }
 
724
                } 
 
725
 
 
726
                /* now solve  u*x = y */
 
727
 
 
728
                for (kb = 0; kb < n; kb++) {
 
729
                    k = n - (kb + 1);
 
730
                    b[k] = b[k]/a[lda*k+k];
 
731
                    t = -b[k];
 
732
                    daxpy(k,t,&a[lda*k+0],1,&b[0],1 );
 
733
                }
 
734
        }
 
735
        else { 
 
736
 
 
737
                /* job = nonzero, solve  trans(a) * x = b
 
738
                   first solve  trans(u)*y = b                  */
 
739
 
 
740
                for (k = 0; k < n; k++) {
 
741
                        t = ddot(k,&a[lda*k+0],1,&b[0],1);
 
742
                        b[k] = (b[k] - t)/a[lda*k+k];
 
743
                }
 
744
 
 
745
                /* now solve trans(l)*x = y     */
 
746
 
 
747
                if (nm1 >= 1) {
 
748
                        for (kb = 1; kb < nm1; kb++) {
 
749
                                k = n - (kb+1);
 
750
                                b[k] = b[k] + ddot(n-(k+1),&a[lda*k+k+1],1,&b[k+1],1);
 
751
                                l = ipvt[k];
 
752
                                if (l != k) {
 
753
                                        t = b[l];
 
754
                                        b[l] = b[k];
 
755
                                        b[k] = t;
 
756
                                }
 
757
                        }
 
758
                }
 
759
        }
 
760
        return;
 
761
}
 
762
 
 
763
/*----------------------*/ 
 
764
 
 
765
void daxpy(int n, REAL da, REAL dx[], int incx, REAL dy[], int incy)
 
766
/*
 
767
     constant times a vector plus a vector.
 
768
     jack dongarra, linpack, 3/11/78.
 
769
*/
 
770
 
 
771
{
 
772
        int i,ix,iy,m,mp1;
 
773
 
 
774
        mp1 = 0;
 
775
        m = 0;
 
776
 
 
777
        if(n <= 0) return;
 
778
        if (da == ZERO) return;
 
779
 
 
780
        if(incx != 1 || incy != 1) {
 
781
 
 
782
                /* code for unequal increments or equal increments
 
783
                   not equal to 1                                       */
 
784
 
 
785
                ix = 0;
 
786
                iy = 0;
 
787
                if(incx < 0) ix = (-n+1)*incx;
 
788
                if(incy < 0)iy = (-n+1)*incy;
 
789
                for (i = 0;i < n; i++) {
 
790
                        dy[iy] = dy[iy] + da*dx[ix];
 
791
                        ix = ix + incx;
 
792
                        iy = iy + incy;
 
793
                     
 
794
                }
 
795
                return;
 
796
        }
 
797
        
 
798
        /* code for both increments equal to 1 */
 
799
        
 
800
 
 
801
#ifdef ROLL
 
802
 
 
803
        for (i = 0;i < n; i++) {
 
804
                dy[i] = dy[i] + da*dx[i];
 
805
        }
 
806
 
 
807
 
 
808
#endif
 
809
 
 
810
#ifdef UNROLL
 
811
 
 
812
        m = n % 4;
 
813
        if ( m != 0) {
 
814
                for (i = 0; i < m; i++) 
 
815
                        dy[i] = dy[i] + da*dx[i];
 
816
                        
 
817
                if (n < 4) return;
 
818
        }
 
819
        for (i = m; i < n; i = i + 4) {
 
820
                dy[i] = dy[i] + da*dx[i];
 
821
                dy[i+1] = dy[i+1] + da*dx[i+1];
 
822
                dy[i+2] = dy[i+2] + da*dx[i+2];
 
823
                dy[i+3] = dy[i+3] + da*dx[i+3];
 
824
                
 
825
        }
 
826
 
 
827
#endif
 
828
return;
 
829
}
 
830
   
 
831
/*----------------------*/ 
 
832
 
 
833
REAL ddot(int n, REAL dx[], int incx, REAL dy[], int incy)
 
834
/*
 
835
     forms the dot product of two vectors.
 
836
     jack dongarra, linpack, 3/11/78.
 
837
*/
 
838
 
 
839
{
 
840
        REAL dtemp;
 
841
        int i,ix,iy,m,mp1;
 
842
 
 
843
        mp1 = 0;
 
844
        m = 0;
 
845
 
 
846
        dtemp = ZERO;
 
847
 
 
848
        if(n <= 0) return(ZERO);
 
849
 
 
850
        if(incx != 1 || incy != 1) {
 
851
 
 
852
                /* code for unequal increments or equal increments
 
853
                   not equal to 1                                       */
 
854
 
 
855
                ix = 0;
 
856
                iy = 0;
 
857
                if (incx < 0) ix = (-n+1)*incx;
 
858
                if (incy < 0) iy = (-n+1)*incy;
 
859
                for (i = 0;i < n; i++) {
 
860
                        dtemp = dtemp + dx[ix]*dy[iy];
 
861
                        ix = ix + incx;
 
862
                        iy = iy + incy;
 
863
                       
 
864
                }
 
865
                return(dtemp);
 
866
        }
 
867
 
 
868
        /* code for both increments equal to 1 */
 
869
 
 
870
 
 
871
#ifdef ROLL
 
872
 
 
873
        for (i=0;i < n; i++)
 
874
                dtemp = dtemp + dx[i]*dy[i];
 
875
               
 
876
        return(dtemp);
 
877
 
 
878
#endif
 
879
 
 
880
#ifdef UNROLL
 
881
 
 
882
 
 
883
        m = n % 5;
 
884
        if (m != 0) {
 
885
                for (i = 0; i < m; i++)
 
886
                        dtemp = dtemp + dx[i]*dy[i];
 
887
                if (n < 5) return(dtemp);
 
888
        }
 
889
        for (i = m; i < n; i = i + 5) {
 
890
                dtemp = dtemp + dx[i]*dy[i] +
 
891
                dx[i+1]*dy[i+1] + dx[i+2]*dy[i+2] +
 
892
                dx[i+3]*dy[i+3] + dx[i+4]*dy[i+4];
 
893
        }
 
894
        return(dtemp);
 
895
 
 
896
#endif
 
897
 
 
898
}
 
899
 
 
900
/*----------------------*/ 
 
901
void dscal(int n, REAL da, REAL dx[], int incx)
 
902
 
 
903
/*     scales a vector by a constant.
 
904
      jack dongarra, linpack, 3/11/78.
 
905
*/
 
906
 
 
907
{
 
908
        int i,m,mp1,nincx;
 
909
 
 
910
        mp1 = 0;
 
911
        m = 0;
 
912
 
 
913
        if(n <= 0)return;
 
914
        if(incx != 1) {
 
915
 
 
916
                /* code for increment not equal to 1 */
 
917
 
 
918
                nincx = n*incx;
 
919
                for (i = 0; i < nincx; i = i + incx)
 
920
                        dx[i] = da*dx[i];
 
921
                        
 
922
                return;
 
923
        }
 
924
 
 
925
        /* code for increment equal to 1 */
 
926
 
 
927
 
 
928
#ifdef ROLL
 
929
 
 
930
        for (i = 0; i < n; i++)
 
931
                dx[i] = da*dx[i];
 
932
                
 
933
 
 
934
#endif
 
935
 
 
936
#ifdef UNROLL
 
937
 
 
938
 
 
939
        m = n % 5;
 
940
        if (m != 0) {
 
941
                for (i = 0; i < m; i++)
 
942
                        dx[i] = da*dx[i];
 
943
                if (n < 5) return;
 
944
        }
 
945
        for (i = m; i < n; i = i + 5){
 
946
                dx[i] = da*dx[i];
 
947
                dx[i+1] = da*dx[i+1];
 
948
                dx[i+2] = da*dx[i+2];
 
949
                dx[i+3] = da*dx[i+3];
 
950
                dx[i+4] = da*dx[i+4];
 
951
        }
 
952
 
 
953
#endif
 
954
 
 
955
}
 
956
 
 
957
/*----------------------*/ 
 
958
int idamax(int n, REAL dx[], int incx)
 
959
 
 
960
/*
 
961
     finds the index of element having max. absolute value.
 
962
     jack dongarra, linpack, 3/11/78.
 
963
*/
 
964
 
 
965
 
 
966
{
 
967
        REAL dmax;
 
968
        int i, ix, itemp;
 
969
 
 
970
        if( n < 1 ) return(-1);
 
971
        if(n ==1 ) return(0);
 
972
        if(incx != 1) {
 
973
 
 
974
                /* code for increment not equal to 1 */
 
975
 
 
976
                ix = 1;
 
977
                dmax = fabs((double)dx[0]);
 
978
                ix = ix + incx;
 
979
                for (i = 1; i < n; i++) {
 
980
                        if(fabs((double)dx[ix]) > dmax)  {
 
981
                                itemp = i;
 
982
                                dmax = fabs((double)dx[ix]);
 
983
                        }
 
984
                        ix = ix + incx;
 
985
                }
 
986
        }
 
987
        else {
 
988
 
 
989
                /* code for increment equal to 1 */
 
990
 
 
991
                itemp = 0;
 
992
                dmax = fabs((double)dx[0]);
 
993
                for (i = 1; i < n; i++) {
 
994
                        if(fabs((double)dx[i]) > dmax) {
 
995
                                itemp = i;
 
996
                                dmax = fabs((double)dx[i]);
 
997
                        }
 
998
                }
 
999
        }
 
1000
        return (itemp);
 
1001
}
 
1002
 
 
1003
/*----------------------*/ 
 
1004
REAL epslon (REAL x)
 
1005
 
 
1006
/*
 
1007
     estimate unit roundoff in quantities of size x.
 
1008
*/
 
1009
 
 
1010
{
 
1011
        REAL a,b,c,eps;
 
1012
/*
 
1013
     this program should function properly on all systems
 
1014
     satisfying the following two assumptions,
 
1015
        1.  the base used in representing dfloating point
 
1016
            numbers is not a power of three.
 
1017
        2.  the quantity  a  in statement 10 is represented to 
 
1018
            the accuracy used in dfloating point variables
 
1019
            that are stored in memory.
 
1020
     the statement number 10 and the go to 10 are intended to
 
1021
     force optimizing compilers to generate code satisfying 
 
1022
     assumption 2.
 
1023
     under these assumptions, it should be true that,
 
1024
            a  is not exactly equal to four-thirds,
 
1025
            b  has a zero for its last bit or digit,
 
1026
            c  is not exactly equal to one,
 
1027
            eps  measures the separation of 1.0 from
 
1028
                 the next larger dfloating point number.
 
1029
     the developers of eispack would appreciate being informed
 
1030
     about any systems where these assumptions do not hold.
 
1031
 
 
1032
     *****************************************************************
 
1033
     this routine is one of the auxiliary routines used by eispack iii
 
1034
     to avoid machine dependencies.
 
1035
     *****************************************************************
 
1036
 
 
1037
     this version dated 4/6/83.
 
1038
*/
 
1039
 
 
1040
        a = 4.0e0/3.0e0;
 
1041
        eps = ZERO;
 
1042
        while (eps == ZERO) {
 
1043
                b = a - ONE;
 
1044
                c = b + b + b;
 
1045
                eps = fabs((double)(c-ONE));
 
1046
        }
 
1047
        return(eps*fabs((double)x));
 
1048
}
 
1049
 
 
1050
/*----------------------*/ 
 
1051
void dmxpy (int n1, REAL y[], int n2, int ldm, REAL x[], REAL m[])
 
1052
 
 
1053
 
 
1054
/* We would like to declare m[][ldm], but c does not allow it.  In this
 
1055
function, references to m[i][j] are written m[ldm*i+j].  */
 
1056
 
 
1057
/*
 
1058
   purpose:
 
1059
     multiply matrix m times vector x and add the result to vector y.
 
1060
 
 
1061
   parameters:
 
1062
 
 
1063
     n1 integer, number of elements in vector y, and number of rows in
 
1064
         matrix m
 
1065
 
 
1066
     y double [n1], vector of length n1 to which is added 
 
1067
         the product m*x
 
1068
 
 
1069
     n2 integer, number of elements in vector x, and number of columns
 
1070
         in matrix m
 
1071
 
 
1072
     ldm integer, leading dimension of array m
 
1073
 
 
1074
     x double [n2], vector of length n2
 
1075
 
 
1076
     m double [ldm][n2], matrix of n1 rows and n2 columns
 
1077
 
 
1078
 ----------------------------------------------------------------------
 
1079
*/
 
1080
{
 
1081
        int j,i,jmin;
 
1082
        /* cleanup odd vector */
 
1083
 
 
1084
        j = n2 % 2;
 
1085
        if (j >= 1) {
 
1086
                j = j - 1;
 
1087
                for (i = 0; i < n1; i++) 
 
1088
                        y[i] = (y[i]) + x[j]*m[ldm*j+i];
 
1089
        } 
 
1090
 
 
1091
        /* cleanup odd group of two vectors */
 
1092
 
 
1093
        j = n2 % 4;
 
1094
        if (j >= 2) {
 
1095
                j = j - 1;
 
1096
                for (i = 0; i < n1; i++)
 
1097
                        y[i] = ( (y[i])
 
1098
                               + x[j-1]*m[ldm*(j-1)+i]) + x[j]*m[ldm*j+i];
 
1099
        } 
 
1100
 
 
1101
        /* cleanup odd group of four vectors */
 
1102
 
 
1103
        j = n2 % 8;
 
1104
        if (j >= 4) {
 
1105
                j = j - 1;
 
1106
                for (i = 0; i < n1; i++)
 
1107
                        y[i] = ((( (y[i])
 
1108
                               + x[j-3]*m[ldm*(j-3)+i]) 
 
1109
                               + x[j-2]*m[ldm*(j-2)+i])
 
1110
                               + x[j-1]*m[ldm*(j-1)+i]) + x[j]*m[ldm*j+i];
 
1111
        } 
 
1112
 
 
1113
        /* cleanup odd group of eight vectors */
 
1114
 
 
1115
        j = n2 % 16;
 
1116
        if (j >= 8) {
 
1117
                j = j - 1;
 
1118
                for (i = 0; i < n1; i++)
 
1119
                        y[i] = ((((((( (y[i])
 
1120
                               + x[j-7]*m[ldm*(j-7)+i]) + x[j-6]*m[ldm*(j-6)+i])
 
1121
                               + x[j-5]*m[ldm*(j-5)+i]) + x[j-4]*m[ldm*(j-4)+i])
 
1122
                               + x[j-3]*m[ldm*(j-3)+i]) + x[j-2]*m[ldm*(j-2)+i])
 
1123
                               + x[j-1]*m[ldm*(j-1)+i]) + x[j]  *m[ldm*j+i];
 
1124
        } 
 
1125
        
 
1126
        /* main loop - groups of sixteen vectors */
 
1127
 
 
1128
        jmin = (n2%16)+16;
 
1129
        for (j = jmin-1; j < n2; j = j + 16) {
 
1130
                for (i = 0; i < n1; i++) 
 
1131
                        y[i] = ((((((((((((((( (y[i])
 
1132
                                + x[j-15]*m[ldm*(j-15)+i]) 
 
1133
                                + x[j-14]*m[ldm*(j-14)+i])
 
1134
                                + x[j-13]*m[ldm*(j-13)+i]) 
 
1135
                                + x[j-12]*m[ldm*(j-12)+i])
 
1136
                                + x[j-11]*m[ldm*(j-11)+i]) 
 
1137
                                + x[j-10]*m[ldm*(j-10)+i])
 
1138
                                + x[j- 9]*m[ldm*(j- 9)+i]) 
 
1139
                                + x[j- 8]*m[ldm*(j- 8)+i])
 
1140
                                + x[j- 7]*m[ldm*(j- 7)+i]) 
 
1141
                                + x[j- 6]*m[ldm*(j- 6)+i])
 
1142
                                + x[j- 5]*m[ldm*(j- 5)+i]) 
 
1143
                                + x[j- 4]*m[ldm*(j- 4)+i])
 
1144
                                + x[j- 3]*m[ldm*(j- 3)+i]) 
 
1145
                                + x[j- 2]*m[ldm*(j- 2)+i])
 
1146
                                + x[j- 1]*m[ldm*(j- 1)+i]) 
 
1147
                                + x[j]   *m[ldm*j+i];
 
1148
        }
 
1149
        return;
 
1150
} 
 
1151
 
 
1152
 
 
1153