← Back to branch summary

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

« back to all changes in this revision

Viewing changes to EXtest/mvttest.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/mvttest.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 transpose, and A always MxN:
 
174
 *   Y is of length N
 
175
 *   X is of length M
 
176
 */
 
177
{
 
178
#ifdef TREAL
 
179
   ATL_INT i, j;
 
180
   TYPE y0;
 
181
 
 
182
   for (j=0; j < N; j++, A += lda)
 
183
   {
 
184
      y0 = Mjoin(PATL,dot)(M, X, incX, A, 1);
 
185
      y0 *= alpha;
 
186
      if (beta == ATL_rone)
 
187
         Y[j] += y0;
 
188
      else if (beta != ATL_rzero)
 
189
         Y[j] = beta*Y[j] + y0;
 
190
      else
 
191
         Y[j] = y0;
 
192
   }
 
193
#else
 
194
   ATL_INT i, j;
 
195
   ATL_CINT lda2 = lda+lda;
 
196
   const TYPE ra=(*alpha), ia=alpha[1], rb=(*beta), ib=beta[1];
 
197
   TYPE ry, iy, rd, id, tmp;
 
198
   TYPE dot[2];
 
199
   void (*mydot)(const int N, const TYPE *X, const int incX,
 
200
                 const TYPE *Y, const int incY, TYPE *dot);
 
201
   mydot = (Conj) ? Mjoin(PATL,dotc_sub) : Mjoin(PATL,dotu_sub);
 
202
 
 
203
   for (j=0; j < N; j++, A += lda2, Y += 2)
 
204
   {
 
205
      mydot(M, X, incX, A, 1, dot);
 
206
/*
 
207
 *    Apply beta to original Y
 
208
 */
 
209
      if (ib == ATL_rzero)
 
210
      {
 
211
         if (rb == ATL_rzero)
 
212
            ry = iy = ATL_rzero;
 
213
         else
 
214
         {
 
215
            ry = *Y * rb;
 
216
            iy = Y[1] * rb;
 
217
         }
 
218
      }
 
219
      else
 
220
      {
 
221
         tmp = ry = *Y;
 
222
         ry = Y[1];
 
223
         ry = ry*rb - iy*ib;
 
224
         iy = tmp*ib + iy*rb;
 
225
      }
 
226
/*
 
227
 *    Apply alpha to dot product
 
228
 */
 
229
      rd = dot[0];
 
230
      id = dot[1];
 
231
      if (ia == ATL_rzero)
 
232
      {
 
233
         ry *= ra;
 
234
         iy *= ra;
 
235
      }
 
236
      else
 
237
      {
 
238
         tmp = ry;
 
239
         ry = ry*ra - iy*ia;
 
240
         iy = tmp*ia + iy*ra;
 
241
      }
 
242
/*
 
243
 *    Store alpha*A*x + beta*y to Y
 
244
 */
 
245
      *Y = ry + rd;
 
246
      Y[1] = iy + id;
 
247
   }
 
248
#endif
 
249
}
 
250
static int CheckAns(int M, int N, TYPE *G, int incG, TYPE *T, int incT)
 
251
{
 
252
   int i, ierr=0;
 
253
   const int NN = N;
 
254
   const double nadds = NN, nmuls = NN+2;
 
255
   #ifdef TREAL
 
256
      const double epsmul = 2*(nadds+nmuls);
 
257
   #else
 
258
      const int incG2 = incG+incG, incT2 = incT+incT;
 
259
      const double epsmul = 2*(nadds+4*nmuls);
 
260
   #endif
 
261
   TYPE maxerr, diff, g, t;
 
262
   maxerr = epsmul * Mjoin(PATL,epsilon)();
 
263
#ifdef TREAL
 
264
   for (i=0; i < NN; i++, G += incG, T += incT)
 
265
   {
 
266
      g = *G;
 
267
      g = Mabs(g);
 
268
      t = *T;
 
269
      t = Mabs(t);
 
270
      diff = g - t;
 
271
      diff = Mabs(diff);
 
272
      if (diff > maxerr)
 
273
      {
 
274
         if (!ierr)
 
275
            ierr = i+1;
 
276
         fprintf(stderr, "Y(%d): Good=%f, Computed=%f, diff=%f\n",
 
277
                 i, *G, *T, diff);
 
278
      }
 
279
   }
 
280
#else
 
281
   for (i=0; i < NN+NN; i += 2, G += incG2, T += incT2)
 
282
   {
 
283
      g = *G;
 
284
      g = Mabs(g);
 
285
      t = *T;
 
286
      t = Mabs(t);
 
287
      diff = g - t;
 
288
      diff = Mabs(diff);
 
289
      if (diff > maxerr)
 
290
      {
 
291
         if (!ierr)
 
292
            ierr = (i>>2)+1;
 
293
         fprintf(stderr, "Yr(%d): Good=%f, Computed=%f, diff=%f\n",
 
294
                 i, *G, *T, diff);
 
295
      }
 
296
      g = G[1];
 
297
      g = Mabs(g);
 
298
      t = T[1];
 
299
      t = Mabs(t);
 
300
      diff = g - t;
 
301
      diff = Mabs(diff);
 
302
      if (diff > maxerr)
 
303
      {
 
304
         if (!ierr)
 
305
            ierr = (i>>2)+1;
 
306
         fprintf(stderr, "Yi(%d): Good=%f, Computed=%f, diff=%f\n",
 
307
                 i, G[1], T[1], diff);
 
308
      }
 
309
   }
 
310
#endif
 
311
   return(ierr);
 
312
}
 
313
#define NX M
 
314
#define NY N
 
315
#ifdef ATL_Conj_
 
316
   #define my_gemv Mjoin(PATL,gemvC)
 
317
#else
 
318
   #define my_gemv Mjoin(PATL,gemvT)
 
319
#endif
 
320
#define ATL_rS2C(sc_) \
 
321
   (((sc_) == ATL_rzero) ? '0' : ( ((sc_) == ATL_rone) ? '1' : 'X'))
 
322
#ifdef TCPLX
 
323
   #define ATL_S2C(sc_) ATL_rS2C(sc_[0]), ATL_rS2C(sc_[1])
 
324
#else
 
325
   #define ATL_S2C(sc_) ATL_rS2C(sc_)
 
326
#endif
 
327
void my_gemv
 
328
   (ATL_CINT M, ATL_CINT N, const SCALAR alpha0, const TYPE *A, ATL_CINT lda,
 
329
    const TYPE *X, ATL_CINT incX, const SCALAR beta0, TYPE *Y, ATL_CINT incY);
 
330
 
 
331
static int RunTest(int M, int N, int lda, int incX, int incY,
 
332
                   SCALAR alpha, SCALAR beta, int II)
 
333
{
 
334
   #ifdef TCPLX
 
335
      TYPE one[2] = {ATL_rone, ATL_rzero};
 
336
   #else
 
337
      TYPE one = ATL_rone;
 
338
   #endif
 
339
   TYPE *A, *A0, *X, *Y, *y;
 
340
   TYPE *Y0;
 
341
   ATL_CINT aincY = Mabs(incY), aincX = Mabs(incX);
 
342
   #ifdef TCPLX
 
343
      char *frm = "%6d %5d %5d %5d %4d %4d %c,%c %c,%c  %4x %4x %4x  %6s\n";
 
344
   #else
 
345
      char *frm = "%6d %5d %5d %5d %4d %4d   %c   %c  %4x %4x %4x  %6s\n";
 
346
   #endif
 
347
   int ierr;
 
348
   if (!II)
 
349
   {
 
350
      printf("\n");
 
351
      printf(
 
352
      "           M     N   lda incY incX alp bet    A    X     Y   PASS?\n");
 
353
      printf(
 
354
      "====== ===== ===== ===== ==== ==== === ===  ==== ==== ====  ======\n");
 
355
   }
 
356
   A = FA_malloc(ATL_MulBySize(lda)*N, FAa, MAa);
 
357
   Y = FA_malloc(ATL_MulBySize(NY)*aincY, FAy, MAy);
 
358
   X = FA_malloc(ATL_MulBySize(NX)*aincX, FAx, MAx);
 
359
   Y0 = FA_malloc(ATL_MulBySize(NY)*aincY, FAy, MAy);
 
360
   ATL_assert(A && Y0 && X && Y);
 
361
   if (aincX != 1)
 
362
      Mjoin(PATLU,set)(NX*aincX SHIFT, 4000000000.0, X, 1);
 
363
   if (aincY != 1)
 
364
   {
 
365
      Mjoin(PATLU,set)(NY*aincY SHIFT, 2000000000.0, Y0, 1);
 
366
      Mjoin(PATLU,set)(NY*aincY SHIFT, 3000000000.0, Y, 1);
 
367
   }
 
368
   Mjoin(PATL,gegen)(1, NY, Y0, aincY, NY*aincY);
 
369
   printf(frm, II, M, N, lda, incY, incX, ATL_S2C(alpha), ATL_S2C(beta),
 
370
          ((size_t)A)&0xFFFF, ((size_t)X)&0xFFFF, ((size_t)Y)&0xFFFF, " START");
 
371
   Mjoin(PATL,gegen)(1, NY, Y, aincY, NY*aincY);
 
372
   Mjoin(PATL,gegen)(1, NX, X, aincX, NY*aincY+127*50+77);
 
373
   Mjoin(PATL,gegen)(M, N, A, lda, N*M+513*7+90);
 
374
   if (incY < 0) Y += (NY-1) * (aincY SHIFT);
 
375
   if (incY < 0) Y0 += (NY-1) * (aincY SHIFT);
 
376
   if (incX < 0) X += (NX-1) * (aincX SHIFT);
 
377
 
 
378
   #ifdef ATL_Conj_
 
379
      dumb_gemv(1, M, N, alpha, A, lda, X, incX, beta, Y0, incY);
 
380
   #else
 
381
      dumb_gemv(0, M, N, alpha, A, lda, X, incX, beta, Y0, incY);
 
382
   #endif
 
383
   my_gemv(M, N, alpha, A, lda, X, incX, beta, Y, incY);
 
384
   ierr = CheckAns(M, N, Y0, incY, Y, incY);
 
385
   FA_free(A, FAa, MAa);
 
386
   if (incY < 0)
 
387
   {
 
388
      Y -= (NY-1) * (aincY SHIFT);
 
389
      Y0 -= (NY-1) * (aincY SHIFT);
 
390
   }
 
391
   FA_free(Y0, FAy, MAy);
 
392
   FA_free(Y, FAy, MAy);
 
393
   if (incX < 0) X -= (NX-1) * (aincX SHIFT);
 
394
   FA_free(X, FAx, MAx);
 
395
   printf(frm, II, M, N, lda, incY, incX, ATL_S2C(alpha), ATL_S2C(beta),
 
396
          ((size_t)A)&0xFFFF, ((size_t) X)&0xFFFF, ((size_t) Y)&0xFFFF,
 
397
          (ierr) ? "FAILED":"PASSED");
 
398
   return(ierr);
 
399
}
 
400
#undef NX
 
401
#undef NY
 
402
 
 
403
int RunTests(int verb, int *Ms, int *Ns, int *incXs, int *incYs, int *ldas,
 
404
             TYPE *alphas, TYPE *betas)
 
405
{
 
406
   int iy, ix, ic, in, im, iax, iay, iaa;
 
407
   ATL_INT m, n, lda, conj, incy;
 
408
   int nerr=0, II=0;
 
409
   int ia, ib, incx;
 
410
   const int NA=(*alphas), NB=(*betas);
 
411
   #ifdef TCPLX
 
412
      TYPE *alpha, *beta;
 
413
   #else
 
414
      TYPE alpha, beta;
 
415
   #endif
 
416
   assert(ldas[0] == Ms[0]);
 
417
   for (in=1; in <= Ns[0]; in++)
 
418
   {
 
419
      n = Ns[in];
 
420
      for (im=1; im <= Ms[0]; im++)
 
421
      {
 
422
         m = Ms[im];
 
423
         lda = ldas[im];
 
424
         for (iy=1; iy <= incYs[0]; iy++)
 
425
         {
 
426
            incy = incYs[iy];
 
427
            for (ix=1; ix <= incXs[0]; ix++)
 
428
            {
 
429
               incx = incXs[ix];
 
430
               for (ia=1; ia <= NA; ia++)
 
431
               {
 
432
                  #ifdef TCPLX
 
433
                     alpha = alphas+ia+ia;
 
434
                  #else
 
435
                     alpha = alphas[ia];
 
436
                  #endif
 
437
                  for (ib=1; ib <= NA; ib++)
 
438
                  {
 
439
                     #ifdef TCPLX
 
440
                        beta = betas+ia+ia;
 
441
                     #else
 
442
                        beta = betas[ia];
 
443
                     #endif
 
444
                     for (iaa=0; iaa < 8; iaa++)
 
445
                     {
 
446
                        FAa = iaa*sizeof(TYPE);
 
447
                        MAa = FAa + sizeof(TYPE);
 
448
                        for (iay=0; iay < 8; iay++)
 
449
                        {
 
450
                           FAy = iay*sizeof(TYPE);
 
451
                           MAy = FAy + sizeof(TYPE);
 
452
                           for (iax=0; iax < 8; iax++)
 
453
                           {
 
454
                              FAx = iax*sizeof(TYPE);
 
455
                              MAx = FAx + sizeof(TYPE);
 
456
                              nerr += RunTest(m, n, lda, incx, incy,
 
457
                                              alpha, beta, II++);
 
458
                           }
 
459
                        }
 
460
                     }
 
461
                  }
 
462
               }
 
463
            }
 
464
         }
 
465
      }
 
466
   }
 
467
   if (nerr)
 
468
      printf("FAILED: %d of %d tests!\n", nerr, II);
 
469
   else
 
470
      printf("PASSED: all %d tests.\n", II);
 
471
   return(nerr);
 
472
}
 
473
 
 
474
void PrintUsage(char *name, int ierr, char *flag)
 
475
{
 
476
   if (ierr > 0)
 
477
      fprintf(stderr, "Bad argument #%d: '%s'\n",
 
478
              ierr, flag ? flag : "Not enough arguments");
 
479
   else if (ierr < 0)
 
480
      fprintf(stderr, "ERROR: %s\n", flag);
 
481
 
 
482
   fprintf(stderr, "USAGE: %s [flags]:\n", name);
 
483
   fprintf(stderr, "   -n <#> <N1> ... <N#>\n");
 
484
   fprintf(stderr, "   -N <Nstart> <Nend> <Ninc>\n");
 
485
   fprintf(stderr, "   -m <#> <M1> ... <M#>\n");
 
486
   fprintf(stderr, "   -M <Mstart> <Mend> <Minc>\n");
 
487
   fprintf(stderr, "   -l <#> <lda1> ... <lda#>\n");
 
488
   fprintf(stderr, "   -g <ldagap> : lda = M + <ldagap> foreach M\n");
 
489
   fprintf(stderr, "   -y <#> <incY1> ... <incY#>\n");
 
490
   fprintf(stderr, "   -x <#> <incX1> ... <incX#>\n");
 
491
   fprintf(stderr, "   -C <#> <conj1> ... <conj#>\n");
 
492
   fprintf(stderr, "   -a <#> <alpha1> ... <alpha#>\n");
 
493
   fprintf(stderr, "   -b <#> <beta1> ... <beta#>\n");
 
494
   fprintf(stderr,
 
495
           "   -v [0,1] : 0 - stop on first error, else keep testing\n");
 
496
   fprintf(stderr,
 
497
"   -F[x,y,a] <#> : if(# > 0) -> force op to be aligned to at least # bytes\n");
 
498
   fprintf(stderr,
 
499
      "                   if(# < 0) -> force op to be aligned to < # bytes.\n");
 
500
 
 
501
   exit(ierr ? ierr : -1);
 
502
}
 
503
 
 
504
 
 
505
/* procedure 1 */
 
506
int *GetIntList1(int ival)
 
507
/*
 
508
 * returns integer array with iarr[0] = 1, iarr[1] = ival
 
509
 */
 
510
{
 
511
   int *iarr;
 
512
   iarr = malloc(2*sizeof(int));
 
513
   ATL_assert(iarr);
 
514
   iarr[0] = 1;
 
515
   iarr[1] = ival;
 
516
   return(iarr);
 
517
}
 
518
 
 
519
#ifdef TYPE
 
520
/* procedure 2 */
 
521
TYPE *GetTypeList1(const SCALAR val)
 
522
/*
 
523
 * Returns a TYPE array with arr[0] = 1.0, arr[1] = val
 
524
 */
 
525
{
 
526
   TYPE *arr;
 
527
   arr = malloc(ATL_MulBySize(2));
 
528
   ATL_assert(arr);
 
529
   arr[0] = 1;
 
530
   #ifdef TCPLX
 
531
      arr[2] = *val;
 
532
      arr[3] = val[1];
 
533
   #else
 
534
      arr[1] = val;
 
535
   #endif
 
536
   return(arr);
 
537
}
 
538
#endif
 
539
 
 
540
/* procedure 3 */
 
541
int *GetIntList2(int ival1, int ival2)
 
542
/*
 
543
 * returns integer array with iarr[0] = 1, iarr[1] = ival1, ival[2] = ival2
 
544
 */
 
545
{
 
546
   int *iarr;
 
547
   iarr = malloc(3*sizeof(int));
 
548
   ATL_assert(iarr);
 
549
   iarr[0] = 1;
 
550
   iarr[1] = ival1;
 
551
   iarr[2] = ival2;
 
552
   return(iarr);
 
553
}
 
554
 
 
555
/* procedure 4 */
 
556
int *DupIntList(int *list)
 
557
/*
 
558
 * Duplicates list of integers, list[0] holds the length, not including 0
 
559
 */
 
560
{
 
561
   int i, n, *ip;
 
562
 
 
563
   assert(list);
 
564
   n = list[0] + 1;
 
565
   ip = malloc(sizeof(int)*n);
 
566
   assert(ip);
 
567
   for (i=0; i < n; i++)
 
568
      ip[i] = list[i];
 
569
   return(ip);
 
570
}
 
571
 
 
572
/* procedure 5 */
 
573
int *GetIntList(int nargs, char **args, int i, int nmul)
 
574
/*
 
575
 * Gets a list of integers, whose length is given by atoi(args[i])*nmul
 
576
 * list is this length+1, since 0'th location gets atoi(args[i])
 
577
 */
 
578
{
 
579
   int n, *iarr, k;
 
580
 
 
581
   if (++i >= nargs)
 
582
      PrintUsage(args[0], i, NULL);
 
583
   n = atoi(args[i]) * nmul;
 
584
   ATL_assert(n > 0);
 
585
   iarr = malloc(sizeof(int)*(n+1));
 
586
   ATL_assert(iarr);
 
587
 
 
588
   iarr[0] = n / nmul;
 
589
   for (k=0; k < n; k++)
 
590
   {
 
591
      if (++i >= nargs)
 
592
         PrintUsage(args[0], i, NULL);
 
593
      iarr[k+1] = atoi(args[i]);
 
594
   }
 
595
   return(iarr);
 
596
}
 
597
 
 
598
#ifdef TYPE
 
599
/* procedure 6 */
 
600
TYPE *GetTypeList(int nargs, char **args, int i, int nmul)
 
601
/*
 
602
 * Gets a list of TYPEs, whose length is given by atoi(args[i])*nmul
 
603
 * list is this length+1, since 0'th location gets atof(args[i])
 
604
 */
 
605
{
 
606
   int n, k;
 
607
   TYPE *arr;
 
608
 
 
609
   if (++i >= nargs)
 
610
      PrintUsage(args[0], i, NULL);
 
611
   n = atoi(args[i]) * nmul;
 
612
   ATL_assert(n > 0);
 
613
   arr = malloc(ATL_MulBySize(n+1));
 
614
   ATL_assert(arr);
 
615
 
 
616
   arr[0] = n / nmul;
 
617
   for (k=0; k < n; k++)
 
618
   {
 
619
      if (++i >= nargs)
 
620
         PrintUsage(args[0], i, NULL);
 
621
      arr[k+(1 SHIFT)] = atof(args[i]);
 
622
   }
 
623
   return(arr);
 
624
}
 
625
#endif
 
626
 
 
627
/* procedure 7 */
 
628
int *IntRange2IntList(int N0, int NN, int incN)
 
629
{
 
630
   int i, n;
 
631
   int *iarr;
 
632
 
 
633
   for (i=N0, n=0; i <= NN; i += incN) n++;
 
634
   iarr = malloc(sizeof(int)*(n+1));
 
635
   ATL_assert(iarr);
 
636
   iarr[0] = n;
 
637
   for (i=N0, n=1 ; i <= NN; i += incN, n++)
 
638
      iarr[n] = i;
 
639
   return(iarr);
 
640
}
 
641
 
 
642
int GetFlags(int nargs, char **args, int **Ms, int **Ns, int **LDAs,
 
643
             int **incYs, int **incXs, TYPE **ALPHAs, TYPE **BETAs)
 
644
{
 
645
   int verb, i, k, *ip;
 
646
   char ch;
 
647
   int ldagap = 8;
 
648
   #ifdef TCPLX
 
649
      const TYPE one[2] = {ATL_rone, ATL_rzero};
 
650
   #else
 
651
      const TYPE one = ATL_rone;
 
652
   #endif
 
653
 
 
654
   *Ns = *Ms = *LDAs = *incYs = *incXs = NULL;
 
655
   *ALPHAs = *BETAs = NULL;
 
656
   verb = 0;
 
657
 
 
658
   for (i=1; i < nargs; i++)
 
659
   {
 
660
      if (args[i][0] != '-')
 
661
         PrintUsage(args[0], i, args[i]);
 
662
      ch = args[i][1];
 
663
      switch(ch)
 
664
      {
 
665
      case 'v':
 
666
        if (++i >= nargs)
 
667
            PrintUsage(args[0], i-1, "out of flags in -g ");
 
668
         verb = atoi(args[i]);
 
669
         break;
 
670
      case 'g':
 
671
        if (++i >= nargs)
 
672
            PrintUsage(args[0], i-1, "out of flags in -g ");
 
673
         ldagap = atoi(args[i]);
 
674
         break;
 
675
      case 'M':
 
676
      case 'N':
 
677
         if (i+3 >= nargs)
 
678
            PrintUsage(args[0], i-1, "out of flags in -N/M ");
 
679
         ip = IntRange2IntList(atoi(args[i+1]),atoi(args[i+2]),atoi(args[i+3]));
 
680
         if (ch == 'M')
 
681
            *Ms = ip;
 
682
         else
 
683
            *Ns = ip;
 
684
         i += 3;
 
685
         break;
 
686
      case 'n':
 
687
      case 'm':
 
688
      case 'l':
 
689
      case 'y':
 
690
      case 'x':
 
691
         ip = GetIntList(nargs, args, i, 1);
 
692
         i += ip[0] + 1;
 
693
         switch(ch)
 
694
         {
 
695
         case 'n':
 
696
            *Ns = ip;
 
697
            break;
 
698
         case 'm':
 
699
            *Ms = ip;
 
700
            break;
 
701
         case 'l':
 
702
            *LDAs = ip;
 
703
            break;
 
704
         case 'y':
 
705
            *incYs = ip;
 
706
            break;
 
707
         case 'x':
 
708
            *incXs = ip;
 
709
            break;
 
710
         }
 
711
         break;
 
712
      case 'a':
 
713
         *ALPHAs = GetTypeList(nargs, args, i, 1 SHIFT);
 
714
         i += (((int)((*ALPHAs)[0]))SHIFT)+1;
 
715
         break;
 
716
      case 'b':
 
717
         *BETAs = GetTypeList(nargs, args, i, 1 SHIFT);
 
718
         i += (((int)((*BETAs)[0]))SHIFT)+1;
 
719
         break;
 
720
      default:
 
721
         PrintUsage(args[0], i, args[i]);
 
722
      }
 
723
   }
 
724
   if (*incXs == NULL)
 
725
      *incXs = GetIntList1(1);
 
726
   if (*incYs == NULL)
 
727
      *incYs = GetIntList1(1);
 
728
   if (*ALPHAs == NULL)
 
729
      *ALPHAs = GetTypeList1(one);
 
730
   if (*BETAs == NULL)
 
731
      *BETAs = GetTypeList1(one);
 
732
   if (*Ms == NULL)
 
733
      *Ms = GetIntList1(977);
 
734
   if (*Ns == NULL)
 
735
      *Ns = GetIntList1(77);
 
736
   if (*LDAs == NULL)
 
737
   {
 
738
      *LDAs = DupIntList(*Ms);
 
739
      for (i=1; i <= (*LDAs)[0]; i++)
 
740
         (*LDAs)[i] += ldagap;
 
741
   }
 
742
   assert((*LDAs)[0] == (*Ms)[0]);
 
743
   return(verb);
 
744
}
 
745
 
 
746
int main(int nargs, char **args)
 
747
{
 
748
   int *Ms, *Ns, *LDAs, *incYs, *incXs, *CONJs;
 
749
   int verb, ierr=0;
 
750
   TYPE *ALPHAs, *BETAs;
 
751
 
 
752
   verb = GetFlags(nargs, args, &Ms, &Ns, &LDAs, &incYs, &incXs,&ALPHAs,&BETAs);
 
753
   ierr = RunTests(verb, Ms, Ns, incXs, incYs, LDAs, ALPHAs, BETAs);
 
754
   free(ALPHAs);
 
755
   free(BETAs);
 
756
   free(incXs);
 
757
   free(incYs);
 
758
   free(Ms);
 
759
   free(Ns);
 
760
   free(LDAs);
 
761
   exit(ierr);
 
762
}