← Back to branch summary

~ubuntu-branches/ubuntu/vivid/atlas/vivid

« back to all changes in this revision

Viewing changes to EXtest/mvntest.c

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-06-11 15:58:16 UTC
  • mfrom: (1.1.3 upstream)
  • mto: (2.2.21 experimental)
  • mto: This revision was merged to the branch mainline in revision 26.
  • Revision ID: package-import@ubuntu.com-20130611155816-b72z8f621tuhbzn0
Tags: upstream-3.10.1
Import upstream version 3.10.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
EXtest/mvntest.c
 
1
/*
 
2
 *             Automatically Tuned Linear Algebra Software v3.10.1
 
3
 * Copyright (C) 2011 R. Clint Whaley
 
4
 *
 
5
 * Redistribution and use in source and binary forms, with or without
 
6
 * modification, are permitted provided that the following conditions
 
7
 * are met:
 
8
 *   1. Redistributions of source code must retain the above copyright
 
9
 *      notice, this list of conditions and the following disclaimer.
 
10
 *   2. Redistributions in binary form must reproduce the above copyright
 
11
 *      notice, this list of conditions, and the following disclaimer in the
 
12
 *      documentation and/or other materials provided with the distribution.
 
13
 *   3. The name of the ATLAS group or the names of its contributers may
 
14
 *      not be used to endorse or promote products derived from this
 
15
 *      software without specific written permission.
 
16
 *
 
17
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 
18
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 
19
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 
20
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
 
21
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 
22
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 
23
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 
24
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 
25
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 
26
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 
27
 * POSSIBILITY OF SUCH DAMAGE.
 
28
 *
 
29
 */
 
30
#include "atlas_misc.h"
 
31
#include "atlas_tst.h"
 
32
#include "atlas_lvl2.h"
 
33
#include "atlas_level1.h"
 
34
#include <ctype.h>
 
35
int FAx=0, MAx=0, FAy=0, MAy=0, FAa=0, MAa=0;
 
36
#include <stdio.h>
 
37
#include <stdlib.h>
 
38
#include <assert.h>
 
39
 
 
40
struct FA_allocs
 
41
{
 
42
   void *mem, *memA;
 
43
   struct FA_allocs *next;
 
44
} *allocQ=NULL;
 
45
 
 
46
struct FA_allocs *NewAlloc(size_t size, struct FA_allocs *next,
 
47
                           int align, int misalign)
 
48
/*
 
49
 * Allocates size allocation that is aligned to [align], but not aligned
 
50
 * to [misalign].  Therefore, misalign > align.  Align must minimally be sizeof
 
51
 * while misalign may be 0 if we don't need to avoid a particular alignment.
 
52
 */
 
53
{
 
54
   void *vp;
 
55
   char *cp;
 
56
   struct FA_allocs *ap;
 
57
   int n, i;
 
58
   const int malign = align >= misalign ? align : misalign;
 
59
 
 
60
   n = size + align + align + malign;
 
61
   i = (n >> 3)<<3;
 
62
   if (n != i)
 
63
      n += n - i;
 
64
   cp = malloc(n + sizeof(struct FA_allocs));
 
65
   assert(cp);
 
66
   ap = (struct FA_allocs *) (cp + n);
 
67
   ap->mem = cp;
 
68
/*
 
69
 * Align to min alignment
 
70
 */
 
71
   ap->memA = align ? (void*) ((((size_t) cp)/align)*align + align) : cp;
 
72
/*
 
73
 * Misalign to misalign
 
74
 * We often need to make sure to unaligned addresses share the same modulo
 
75
 * so that they have the *same* degree of misalignment (so that their alignment
 
76
 * can be fixed by simple peeling), and so in this case force the address
 
77
 * modulo the misalign to be the exact align value.
 
78
 */
 
79
   if (misalign)
 
80
      ap->memA = (void*)((((size_t)ap->memA)/malign)*malign + malign + align);
 
81
   ap->next = next;
 
82
   return(ap);
 
83
}
 
84
 
 
85
/*
 
86
 * no-align malloc free retaining system default behavior
 
87
 */
 
88
void *NA_malloc(size_t size)
 
89
{
 
90
   return(malloc(size));
 
91
}
 
92
void *NA_calloc(size_t n, size_t size)
 
93
{
 
94
   return(calloc(n, size));
 
95
}
 
96
void NA_free(void *ptr)
 
97
{
 
98
   free(ptr);
 
99
}
 
100
 
 
101
 
 
102
/*
 
103
 * malloc/free pair that aligns data to align, but not to misalign
 
104
 */
 
105
void *FA_malloc(size_t size, int align, int misalign)
 
106
{
 
107
   if ((!misalign && align <= 8) || !size)
 
108
      return(malloc(size));
 
109
   else
 
110
   {
 
111
      allocQ = NewAlloc(size, allocQ, align, misalign);
 
112
      return(allocQ->memA);
 
113
   }
 
114
}
 
115
void *FA_calloc(size_t n, size_t size, int align, int misalign)
 
116
{
 
117
   char *cp;
 
118
   int *ip;
 
119
   double *dp;
 
120
   size_t i;
 
121
   size_t tsize;
 
122
   tsize = n * size;
 
123
   cp = FA_malloc(tsize, align, misalign);
 
124
   if (size == sizeof(int))
 
125
      for (ip=(int*)cp,i=0; i < n; i++)
 
126
        ip[i] = 0;
 
127
   else if (size == sizeof(double))
 
128
      for (dp=(double*)cp,i=0; i < n; i++)
 
129
        dp[i] = 0.0;
 
130
   else
 
131
      for (i=0; i < tsize; i++)
 
132
        cp[i] = 0;
 
133
   return(cp);
 
134
}
 
135
 
 
136
void FA_free(void *ptr, int align, int misalign)
 
137
/*
 
138
 * Part of malloc/free pair that aligns data to FALIGN
 
139
 */
 
140
{
 
141
   struct FA_allocs *ap, *prev;
 
142
   if (ptr)
 
143
   {
 
144
      if ((!misalign && align <= 8))
 
145
         free(ptr);
 
146
      else
 
147
      {
 
148
         for (ap=allocQ; ap && ap->memA != ptr; ap = ap->next) prev = ap;
 
149
         if (!ap)
 
150
         {
 
151
            fprintf(stderr, "Couldn't find mem=%ld\nmemQ=\n", (size_t)ptr);
 
152
            for (ap=allocQ; ap; ap = ap->next)
 
153
               fprintf(stderr, "   %ld, %ld\n", (size_t)ap->memA,
 
154
                       (size_t)ap->mem);
 
155
         }
 
156
         assert(ap);
 
157
         if (ap == allocQ)
 
158
            allocQ = allocQ->next;
 
159
         else
 
160
            prev->next = ap->next;
 
161
         free(ap->mem);
 
162
      }
 
163
   }
 
164
}
 
165
 
 
166
static void dumb_gemv
 
167
   (int Conj, ATL_CINT M, ATL_CINT N, const SCALAR alpha,
 
168
    const TYPE *A, ATL_CINT lda, const TYPE *X, ATL_CINT incX,
 
169
    const SCALAR beta, TYPE *Y, ATL_CINT incY)
 
170
/*
 
171
 * If (conj) y = alpha*conj(A^T)*x + beta*y
 
172
 * else  y = alpha*(A^T)*x + beta*y
 
173
 * Since this is no-transpose, and A always MxN:
 
174
 *   Y is of length M
 
175
 *   X is of length N
 
176
 */
 
177
{
 
178
   ATL_INT j;
 
179
   #ifdef TCPLX
 
180
      const TYPE ra=(*alpha), ia=alpha[1];
 
181
      TYPE rx, ix, xa[2];
 
182
      const int lda2 = lda+lda, incX2 = incX+incX;
 
183
   #else
 
184
      TYPE xa;
 
185
      const int lda2 = lda, incX2 = incX;
 
186
   #endif
 
187
/*
 
188
 * Apply beta to Y
 
189
 */
 
190
   if (!SCALAR_IS_ONE(beta))
 
191
   {
 
192
      if (SCALAR_IS_ZERO(beta))
 
193
         Mjoin(PATL,zero)(M, Y, incY);
 
194
      else
 
195
         Mjoin(PATL,scal)(M, beta, Y, incY);
 
196
   }
 
197
   if (SCALAR_IS_ZERO(alpha) || N < 1)
 
198
      return;
 
199
   for (j=0; j < N; j++, A += lda2, X += incX2)
 
200
   {
 
201
      #ifdef TCPLX
 
202
         rx = *X;
 
203
         ix = X[1];
 
204
         xa[0] = rx*ra - ix*ia;
 
205
         xa[1] = rx*ia + ix*ra;
 
206
      #else
 
207
        xa = alpha * *X;
 
208
      #endif
 
209
      Mjoin(PATL,axpy)(M, xa, A, 1, Y, incY);
 
210
   }
 
211
}
 
212
static int CheckAns(int M, int N, TYPE *G, int incG, TYPE *T, int incT)
 
213
{
 
214
   int i, ierr=0;
 
215
   const int NN = M;
 
216
   const double nadds = NN, nmuls = NN+2;
 
217
   #ifdef TREAL
 
218
      const double epsmul = 2*(nadds+nmuls);
 
219
   #else
 
220
      const int incG2 = incG+incG, incT2 = incT+incT;
 
221
      const double epsmul = 2*(nadds+4*nmuls);
 
222
   #endif
 
223
   TYPE maxerr, diff, g, t;
 
224
   maxerr = epsmul * Mjoin(PATL,epsilon)();
 
225
#ifdef TREAL
 
226
   for (i=0; i < NN; i++, G += incG, T += incT)
 
227
   {
 
228
      g = *G;
 
229
      g = Mabs(g);
 
230
      t = *T;
 
231
      t = Mabs(t);
 
232
      diff = g - t;
 
233
      diff = Mabs(diff);
 
234
      if (diff > maxerr)
 
235
      {
 
236
         if (!ierr)
 
237
            ierr = i+1;
 
238
         fprintf(stderr, "Y(%d): Good=%f, Computed=%f, diff=%f\n",
 
239
                 i, *G, *T, diff);
 
240
      }
 
241
   }
 
242
#else
 
243
   for (i=0; i < NN+NN; i += 2, G += incG2, T += incT2)
 
244
   {
 
245
      g = *G;
 
246
      g = Mabs(g);
 
247
      t = *T;
 
248
      t = Mabs(t);
 
249
      diff = g - t;
 
250
      diff = Mabs(diff);
 
251
      if (diff > maxerr)
 
252
      {
 
253
         if (!ierr)
 
254
            ierr = (i>>2)+1;
 
255
         fprintf(stderr, "Yr(%d): Good=%f, Computed=%f, diff=%f\n",
 
256
                 i, *G, *T, diff);
 
257
      }
 
258
      g = G[1];
 
259
      g = Mabs(g);
 
260
      t = T[1];
 
261
      t = Mabs(t);
 
262
      diff = g - t;
 
263
      diff = Mabs(diff);
 
264
      if (diff > maxerr)
 
265
      {
 
266
         if (!ierr)
 
267
            ierr = (i>>2)+1;
 
268
         fprintf(stderr, "Yi(%d): Good=%f, Computed=%f, diff=%f\n",
 
269
                 i, G[1], T[1], diff);
 
270
      }
 
271
   }
 
272
#endif
 
273
   return(ierr);
 
274
}
 
275
#define NX N
 
276
#define NY M
 
277
#define my_gemv Mjoin(PATL,gemvN)
 
278
#define ATL_rS2C(sc_) \
 
279
   (((sc_) == ATL_rzero) ? '0' : ( ((sc_) == ATL_rone) ? '1' : 'X'))
 
280
#ifdef TCPLX
 
281
   #define ATL_S2C(sc_) ATL_rS2C(sc_[0]), ATL_rS2C(sc_[1])
 
282
#else
 
283
   #define ATL_S2C(sc_) ATL_rS2C(sc_)
 
284
#endif
 
285
void my_gemv
 
286
   (ATL_CINT M, ATL_CINT N, const SCALAR alpha0, const TYPE *A, ATL_CINT lda,
 
287
    const TYPE *X, ATL_CINT incX, const SCALAR beta0, TYPE *Y, ATL_CINT incY);
 
288
 
 
289
static int RunTest(int M, int N, int lda, int incX, int incY,
 
290
                   SCALAR alpha, SCALAR beta, int II)
 
291
{
 
292
   #ifdef TCPLX
 
293
      TYPE one[2] = {ATL_rone, ATL_rzero};
 
294
   #else
 
295
      TYPE one = ATL_rone;
 
296
   #endif
 
297
   TYPE *A, *A0, *X, *Y, *y;
 
298
   TYPE *Y0;
 
299
   ATL_CINT aincY = Mabs(incY), aincX = Mabs(incX);
 
300
   #ifdef TCPLX
 
301
      char *frm = "%6d %5d %5d %5d %4d %4d %c,%c %c,%c  %4x %4x %4x  %6s\n";
 
302
   #else
 
303
      char *frm = "%6d %5d %5d %5d %4d %4d   %c   %c  %4x %4x %4x  %6s\n";
 
304
   #endif
 
305
   int ierr;
 
306
   if (!II)
 
307
   {
 
308
      printf("\n");
 
309
      printf(
 
310
      "           M     N   lda incY incX alp bet    A    X     Y   PASS?\n");
 
311
      printf(
 
312
      "====== ===== ===== ===== ==== ==== === ===  ==== ==== ====  ======\n");
 
313
   }
 
314
   A = FA_malloc(ATL_MulBySize(lda)*N, FAa, MAa);
 
315
   Y = FA_malloc(ATL_MulBySize(NY)*aincY, FAy, MAy);
 
316
   X = FA_malloc(ATL_MulBySize(NX)*aincX, FAx, MAx);
 
317
   Y0 = FA_malloc(ATL_MulBySize(NY)*aincY, FAy, MAy);
 
318
   ATL_assert(A && Y0 && X && Y);
 
319
   if (aincX != 1)
 
320
      Mjoin(PATLU,set)(NX*aincX SHIFT, 4000000000.0, X, 1);
 
321
   if (aincY != 1)
 
322
   {
 
323
      Mjoin(PATLU,set)(NY*aincY SHIFT, 2000000000.0, Y0, 1);
 
324
      Mjoin(PATLU,set)(NY*aincY SHIFT, 3000000000.0, Y, 1);
 
325
   }
 
326
   Mjoin(PATL,gegen)(1, NY, Y0, aincY, NY*aincY);
 
327
   printf(frm, II, M, N, lda, incY, incX, ATL_S2C(alpha), ATL_S2C(beta),
 
328
          ((size_t)A)&0xFFFF, ((size_t)X)&0xFFFF, ((size_t)Y)&0xFFFF, " START");
 
329
   Mjoin(PATL,gegen)(1, NY, Y, aincY, NY*aincY);
 
330
   Mjoin(PATL,gegen)(1, NX, X, aincX, NY*aincY+127*50+77);
 
331
   Mjoin(PATL,gegen)(M, N, A, lda, N*M+513*7+90);
 
332
   if (incY < 0) Y += (NY-1) * (aincY SHIFT);
 
333
   if (incY < 0) Y0 += (NY-1) * (aincY SHIFT);
 
334
   if (incX < 0) X += (NX-1) * (aincX SHIFT);
 
335
 
 
336
   #ifdef ATL_Conj_
 
337
      dumb_gemv(1, M, N, alpha, A, lda, X, incX, beta, Y0, incY);
 
338
   #else
 
339
      dumb_gemv(0, M, N, alpha, A, lda, X, incX, beta, Y0, incY);
 
340
   #endif
 
341
   my_gemv(M, N, alpha, A, lda, X, incX, beta, Y, incY);
 
342
   ierr = CheckAns(M, N, Y0, incY, Y, incY);
 
343
   FA_free(A, FAa, MAa);
 
344
   if (incY < 0)
 
345
   {
 
346
      Y -= (NY-1) * (aincY SHIFT);
 
347
      Y0 -= (NY-1) * (aincY SHIFT);
 
348
   }
 
349
   FA_free(Y0, FAy, MAy);
 
350
   FA_free(Y, FAy, MAy);
 
351
   if (incX < 0) X -= (NX-1) * (aincX SHIFT);
 
352
   FA_free(X, FAx, MAx);
 
353
   printf(frm, II, M, N, lda, incY, incX, ATL_S2C(alpha), ATL_S2C(beta),
 
354
          ((size_t)A)&0xFFFF, ((size_t) X)&0xFFFF, ((size_t) Y)&0xFFFF,
 
355
          (ierr) ? "FAILED":"PASSED");
 
356
   return(ierr);
 
357
}
 
358
#undef NX
 
359
#undef NY
 
360
 
 
361
int RunTests(int verb, int *Ms, int *Ns, int *incXs, int *incYs, int *ldas,
 
362
             TYPE *alphas, TYPE *betas)
 
363
{
 
364
   int iy, ix, ic, in, im, iax, iay, iaa;
 
365
   ATL_INT m, n, lda, conj, incy;
 
366
   int nerr=0, II=0;
 
367
   int ia, ib, incx;
 
368
   const int NA=(*alphas), NB=(*betas);
 
369
   #ifdef TCPLX
 
370
      TYPE *alpha, *beta;
 
371
   #else
 
372
      TYPE alpha, beta;
 
373
   #endif
 
374
   assert(ldas[0] == Ms[0]);
 
375
   for (in=1; in <= Ns[0]; in++)
 
376
   {
 
377
      n = Ns[in];
 
378
      for (im=1; im <= Ms[0]; im++)
 
379
      {
 
380
         m = Ms[im];
 
381
         lda = ldas[im];
 
382
         for (iy=1; iy <= incYs[0]; iy++)
 
383
         {
 
384
            incy = incYs[iy];
 
385
            for (ix=1; ix <= incXs[0]; ix++)
 
386
            {
 
387
               incx = incXs[ix];
 
388
               for (ia=1; ia <= NA; ia++)
 
389
               {
 
390
                  #ifdef TCPLX
 
391
                     alpha = alphas+ia+ia;
 
392
                  #else
 
393
                     alpha = alphas[ia];
 
394
                  #endif
 
395
                  for (ib=1; ib <= NA; ib++)
 
396
                  {
 
397
                     #ifdef TCPLX
 
398
                        beta = betas+ia+ia;
 
399
                     #else
 
400
                        beta = betas[ia];
 
401
                     #endif
 
402
                     for (iaa=0; iaa < 8; iaa++)
 
403
                     {
 
404
                        FAa = iaa*sizeof(TYPE);
 
405
                        MAa = FAa + sizeof(TYPE);
 
406
                        for (iay=0; iay < 8; iay++)
 
407
                        {
 
408
                           FAy = iay*sizeof(TYPE);
 
409
                           MAy = FAy + sizeof(TYPE);
 
410
                           for (iax=0; iax < 8; iax++)
 
411
                           {
 
412
                              FAx = iax*sizeof(TYPE);
 
413
                              MAx = FAx + sizeof(TYPE);
 
414
                              nerr += RunTest(m, n, lda, incx, incy,
 
415
                                              alpha, beta, II++);
 
416
                           }
 
417
                        }
 
418
                     }
 
419
                  }
 
420
               }
 
421
            }
 
422
         }
 
423
      }
 
424
   }
 
425
   if (nerr)
 
426
      printf("FAILED: %d of %d tests!\n", nerr, II);
 
427
   else
 
428
      printf("PASSED: all %d tests.\n", II);
 
429
   return(nerr);
 
430
}
 
431
 
 
432
void PrintUsage(char *name, int ierr, char *flag)
 
433
{
 
434
   if (ierr > 0)
 
435
      fprintf(stderr, "Bad argument #%d: '%s'\n",
 
436
              ierr, flag ? flag : "Not enough arguments");
 
437
   else if (ierr < 0)
 
438
      fprintf(stderr, "ERROR: %s\n", flag);
 
439
 
 
440
   fprintf(stderr, "USAGE: %s [flags]:\n", name);
 
441
   fprintf(stderr, "   -n <#> <N1> ... <N#>\n");
 
442
   fprintf(stderr, "   -N <Nstart> <Nend> <Ninc>\n");
 
443
   fprintf(stderr, "   -m <#> <M1> ... <M#>\n");
 
444
   fprintf(stderr, "   -M <Mstart> <Mend> <Minc>\n");
 
445
   fprintf(stderr, "   -l <#> <lda1> ... <lda#>\n");
 
446
   fprintf(stderr, "   -g <ldagap> : lda = M + <ldagap> foreach M\n");
 
447
   fprintf(stderr, "   -y <#> <incY1> ... <incY#>\n");
 
448
   fprintf(stderr, "   -x <#> <incX1> ... <incX#>\n");
 
449
   fprintf(stderr, "   -C <#> <conj1> ... <conj#>\n");
 
450
   fprintf(stderr, "   -a <#> <alpha1> ... <alpha#>\n");
 
451
   fprintf(stderr, "   -b <#> <beta1> ... <beta#>\n");
 
452
   fprintf(stderr,
 
453
           "   -v [0,1] : 0 - stop on first error, else keep testing\n");
 
454
   fprintf(stderr,
 
455
"   -F[x,y,a] <#> : if(# > 0) -> force op to be aligned to at least # bytes\n");
 
456
   fprintf(stderr,
 
457
      "                   if(# < 0) -> force op to be aligned to < # bytes.\n");
 
458
 
 
459
   exit(ierr ? ierr : -1);
 
460
}
 
461
 
 
462
 
 
463
/* procedure 1 */
 
464
int *GetIntList1(int ival)
 
465
/*
 
466
 * returns integer array with iarr[0] = 1, iarr[1] = ival
 
467
 */
 
468
{
 
469
   int *iarr;
 
470
   iarr = malloc(2*sizeof(int));
 
471
   ATL_assert(iarr);
 
472
   iarr[0] = 1;
 
473
   iarr[1] = ival;
 
474
   return(iarr);
 
475
}
 
476
 
 
477
#ifdef TYPE
 
478
/* procedure 2 */
 
479
TYPE *GetTypeList1(const SCALAR val)
 
480
/*
 
481
 * Returns a TYPE array with arr[0] = 1.0, arr[1] = val
 
482
 */
 
483
{
 
484
   TYPE *arr;
 
485
   arr = malloc(ATL_MulBySize(2));
 
486
   ATL_assert(arr);
 
487
   arr[0] = 1;
 
488
   #ifdef TCPLX
 
489
      arr[2] = *val;
 
490
      arr[3] = val[1];
 
491
   #else
 
492
      arr[1] = val;
 
493
   #endif
 
494
   return(arr);
 
495
}
 
496
#endif
 
497
 
 
498
/* procedure 3 */
 
499
int *GetIntList2(int ival1, int ival2)
 
500
/*
 
501
 * returns integer array with iarr[0] = 1, iarr[1] = ival1, ival[2] = ival2
 
502
 */
 
503
{
 
504
   int *iarr;
 
505
   iarr = malloc(3*sizeof(int));
 
506
   ATL_assert(iarr);
 
507
   iarr[0] = 1;
 
508
   iarr[1] = ival1;
 
509
   iarr[2] = ival2;
 
510
   return(iarr);
 
511
}
 
512
 
 
513
/* procedure 4 */
 
514
int *DupIntList(int *list)
 
515
/*
 
516
 * Duplicates list of integers, list[0] holds the length, not including 0
 
517
 */
 
518
{
 
519
   int i, n, *ip;
 
520
 
 
521
   assert(list);
 
522
   n = list[0] + 1;
 
523
   ip = malloc(sizeof(int)*n);
 
524
   assert(ip);
 
525
   for (i=0; i < n; i++)
 
526
      ip[i] = list[i];
 
527
   return(ip);
 
528
}
 
529
 
 
530
/* procedure 5 */
 
531
int *GetIntList(int nargs, char **args, int i, int nmul)
 
532
/*
 
533
 * Gets a list of integers, whose length is given by atoi(args[i])*nmul
 
534
 * list is this length+1, since 0'th location gets atoi(args[i])
 
535
 */
 
536
{
 
537
   int n, *iarr, k;
 
538
 
 
539
   if (++i >= nargs)
 
540
      PrintUsage(args[0], i, NULL);
 
541
   n = atoi(args[i]) * nmul;
 
542
   ATL_assert(n > 0);
 
543
   iarr = malloc(sizeof(int)*(n+1));
 
544
   ATL_assert(iarr);
 
545
 
 
546
   iarr[0] = n / nmul;
 
547
   for (k=0; k < n; k++)
 
548
   {
 
549
      if (++i >= nargs)
 
550
         PrintUsage(args[0], i, NULL);
 
551
      iarr[k+1] = atoi(args[i]);
 
552
   }
 
553
   return(iarr);
 
554
}
 
555
 
 
556
#ifdef TYPE
 
557
/* procedure 6 */
 
558
TYPE *GetTypeList(int nargs, char **args, int i, int nmul)
 
559
/*
 
560
 * Gets a list of TYPEs, whose length is given by atoi(args[i])*nmul
 
561
 * list is this length+1, since 0'th location gets atof(args[i])
 
562
 */
 
563
{
 
564
   int n, k;
 
565
   TYPE *arr;
 
566
 
 
567
   if (++i >= nargs)
 
568
      PrintUsage(args[0], i, NULL);
 
569
   n = atoi(args[i]) * nmul;
 
570
   ATL_assert(n > 0);
 
571
   arr = malloc(ATL_MulBySize(n+1));
 
572
   ATL_assert(arr);
 
573
 
 
574
   arr[0] = n / nmul;
 
575
   for (k=0; k < n; k++)
 
576
   {
 
577
      if (++i >= nargs)
 
578
         PrintUsage(args[0], i, NULL);
 
579
      arr[k+(1 SHIFT)] = atof(args[i]);
 
580
   }
 
581
   return(arr);
 
582
}
 
583
#endif
 
584
 
 
585
/* procedure 7 */
 
586
int *IntRange2IntList(int N0, int NN, int incN)
 
587
{
 
588
   int i, n;
 
589
   int *iarr;
 
590
 
 
591
   for (i=N0, n=0; i <= NN; i += incN) n++;
 
592
   iarr = malloc(sizeof(int)*(n+1));
 
593
   ATL_assert(iarr);
 
594
   iarr[0] = n;
 
595
   for (i=N0, n=1 ; i <= NN; i += incN, n++)
 
596
      iarr[n] = i;
 
597
   return(iarr);
 
598
}
 
599
 
 
600
int GetFlags(int nargs, char **args, int **Ms, int **Ns, int **LDAs,
 
601
             int **incYs, int **incXs, TYPE **ALPHAs, TYPE **BETAs)
 
602
{
 
603
   int verb, i, k, *ip;
 
604
   char ch;
 
605
   int ldagap = 8;
 
606
   #ifdef TCPLX
 
607
      const TYPE one[2] = {ATL_rone, ATL_rzero};
 
608
   #else
 
609
      const TYPE one = ATL_rone;
 
610
   #endif
 
611
 
 
612
   *Ns = *Ms = *LDAs = *incYs = *incXs = NULL;
 
613
   *ALPHAs = *BETAs = NULL;
 
614
   verb = 0;
 
615
 
 
616
   for (i=1; i < nargs; i++)
 
617
   {
 
618
      if (args[i][0] != '-')
 
619
         PrintUsage(args[0], i, args[i]);
 
620
      ch = args[i][1];
 
621
      switch(ch)
 
622
      {
 
623
      case 'v':
 
624
        if (++i >= nargs)
 
625
            PrintUsage(args[0], i-1, "out of flags in -g ");
 
626
         verb = atoi(args[i]);
 
627
         break;
 
628
      case 'g':
 
629
        if (++i >= nargs)
 
630
            PrintUsage(args[0], i-1, "out of flags in -g ");
 
631
         ldagap = atoi(args[i]);
 
632
         break;
 
633
      case 'M':
 
634
      case 'N':
 
635
         if (i+3 >= nargs)
 
636
            PrintUsage(args[0], i-1, "out of flags in -N/M ");
 
637
         ip = IntRange2IntList(atoi(args[i+1]),atoi(args[i+2]),atoi(args[i+3]));
 
638
         if (ch == 'M')
 
639
            *Ms = ip;
 
640
         else
 
641
            *Ns = ip;
 
642
         i += 3;
 
643
         break;
 
644
      case 'n':
 
645
      case 'm':
 
646
      case 'l':
 
647
      case 'y':
 
648
      case 'x':
 
649
         ip = GetIntList(nargs, args, i, 1);
 
650
         i += ip[0] + 1;
 
651
         switch(ch)
 
652
         {
 
653
         case 'n':
 
654
            *Ns = ip;
 
655
            break;
 
656
         case 'm':
 
657
            *Ms = ip;
 
658
            break;
 
659
         case 'l':
 
660
            *LDAs = ip;
 
661
            break;
 
662
         case 'y':
 
663
            *incYs = ip;
 
664
            break;
 
665
         case 'x':
 
666
            *incXs = ip;
 
667
            break;
 
668
         }
 
669
         break;
 
670
      case 'a':
 
671
         *ALPHAs = GetTypeList(nargs, args, i, 1 SHIFT);
 
672
         i += (((int)((*ALPHAs)[0]))SHIFT)+1;
 
673
         break;
 
674
      case 'b':
 
675
         *BETAs = GetTypeList(nargs, args, i, 1 SHIFT);
 
676
         i += (((int)((*BETAs)[0]))SHIFT)+1;
 
677
         break;
 
678
      default:
 
679
         PrintUsage(args[0], i, args[i]);
 
680
      }
 
681
   }
 
682
   if (*incXs == NULL)
 
683
      *incXs = GetIntList1(1);
 
684
   if (*incYs == NULL)
 
685
      *incYs = GetIntList1(1);
 
686
   if (*ALPHAs == NULL)
 
687
      *ALPHAs = GetTypeList1(one);
 
688
   if (*BETAs == NULL)
 
689
      *BETAs = GetTypeList1(one);
 
690
   if (*Ms == NULL)
 
691
      *Ms = GetIntList1(977);
 
692
   if (*Ns == NULL)
 
693
      *Ns = GetIntList1(77);
 
694
   if (*LDAs == NULL)
 
695
   {
 
696
      *LDAs = DupIntList(*Ms);
 
697
      for (i=1; i <= (*LDAs)[0]; i++)
 
698
         (*LDAs)[i] += ldagap;
 
699
   }
 
700
   assert((*LDAs)[0] == (*Ms)[0]);
 
701
   return(verb);
 
702
}
 
703
 
 
704
int main(int nargs, char **args)
 
705
{
 
706
   int *Ms, *Ns, *LDAs, *incYs, *incXs, *CONJs;
 
707
   int verb, ierr=0;
 
708
   TYPE *ALPHAs, *BETAs;
 
709
 
 
710
   verb = GetFlags(nargs, args, &Ms, &Ns, &LDAs, &incYs, &incXs,&ALPHAs,&BETAs);
 
711
   ierr = RunTests(verb, Ms, Ns, incXs, incYs, LDAs, ALPHAs, BETAs);
 
712
   free(ALPHAs);
 
713
   free(BETAs);
 
714
   free(incXs);
 
715
   free(incYs);
 
716
   free(Ms);
 
717
   free(Ns);
 
718
   free(LDAs);
 
719
   exit(ierr);
 
720
}