memtort.c

来自「C语言版本的矩阵库」· C语言 代码 · 共 761 行 · 第 1/2 页

/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**			     Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/


/* 
  Tests for mem_info.c functions
  */

static char rcsid[] = "$Id: $";

#include        <stdio.h>
#include        <math.h>
#include        "matrix2.h"
#include 	"sparse2.h"
#include  	"zmatrix2.h"


#define errmesg(mesg)   printf("Error: %s error: line %d\n",mesg,__LINE__)
#define notice(mesg)    printf("# Testing %s...\n",mesg)


/*  new types list */

extern MEM_CONNECT mem_connect[MEM_CONNECT_MAX_LISTS];

/* the number of a new list */
#define FOO_LIST 1

/* numbers of types */
#define TYPE_FOO_1    1
#define TYPE_FOO_2    2

typedef struct {
   int dim;
   int fix_dim;
   Real (*a)[10];
} FOO_1;

typedef struct {
  int dim;
  int fix_dim;
  Real (*a)[2];
} FOO_2;



FOO_1 *foo_1_get(dim)
int dim;
{
   FOO_1 *f;
   
   if ((f = (FOO_1 *)malloc(sizeof(FOO_1))) == NULL)
     error(E_MEM,"foo_1_get");
   else if (mem_info_is_on()) {
      mem_bytes_list(TYPE_FOO_1,0,sizeof(FOO_1),FOO_LIST);
      mem_numvar_list(TYPE_FOO_1,1,FOO_LIST);
   }
   
   f->dim = dim;
   f->fix_dim = 10;
   if ((f->a = (Real (*)[10])malloc(dim*sizeof(Real [10]))) == NULL)
      error(E_MEM,"foo_1_get");
   else if (mem_info_is_on())
     mem_bytes_list(TYPE_FOO_1,0,dim*sizeof(Real [10]),FOO_LIST); 

   return f;
}


FOO_2 *foo_2_get(dim)
int dim;
{
   FOO_2 *f;
   
   if ((f = (FOO_2 *)malloc(sizeof(FOO_2))) == NULL)
     error(E_MEM,"foo_2_get");
   else if (mem_info_is_on()) {
      mem_bytes_list(TYPE_FOO_2,0,sizeof(FOO_2),FOO_LIST);
      mem_numvar_list(TYPE_FOO_2,1,FOO_LIST);
   }

   f->dim = dim;
   f->fix_dim = 2;
   if ((f->a = (Real (*)[2])malloc(dim*sizeof(Real [2]))) == NULL)
      error(E_MEM,"foo_2_get");
   else if (mem_info_is_on())
     mem_bytes_list(TYPE_FOO_2,0,dim*sizeof(Real [2]),FOO_LIST); 

   return f;
}



int foo_1_free(f)
FOO_1 *f;
{
   if ( f != NULL) {
      if (mem_info_is_on()) {
	 mem_bytes_list(TYPE_FOO_1,sizeof(FOO_1)+
			f->dim*sizeof(Real [10]),0,FOO_LIST);
	 mem_numvar_list(TYPE_FOO_1,-1,FOO_LIST);
      }

      free(f->a);
      free(f);
   }
   return 0;
}

int foo_2_free(f)
FOO_2 *f;
{
   if ( f != NULL) {
      if (mem_info_is_on()) {
	 mem_bytes_list(TYPE_FOO_2,sizeof(FOO_2)+
			f->dim*sizeof(Real [2]),0,FOO_LIST);
	 mem_numvar_list(TYPE_FOO_2,-1,FOO_LIST);
      }

      free(f->a);
      free(f);
   }
   return 0;
}




char *foo_type_name[] = {
   "nothing",
   "FOO_1",
   "FOO_2"
};


#define FOO_NUM_TYPES  (sizeof(foo_type_name)/sizeof(*foo_type_name))


int (*foo_free_func[FOO_NUM_TYPES])() = {
   NULL, 
   foo_1_free, 
   foo_2_free
  };



static MEM_ARRAY foo_info_sum[FOO_NUM_TYPES];



  /* px_rand -- generates sort-of random permutation */
PERM    *px_rand(pi)
PERM    *pi;
{
   int         i, j, k;
   
   if ( ! pi )
     error(E_NULL,"px_rand");
   
   for ( i = 0; i < 3*pi->size; i++ )
   {
      j = (rand() >> 8) % pi->size;
      k = (rand() >> 8) % pi->size;
      px_transp(pi,j,k);
   }
   
   return pi;
}

#ifdef SPARSE
SPMAT  *gen_non_symm(m,n)
int     m, n;
{
    SPMAT      *A;
    static      PERM    *px = PNULL;
    int         i, j, k, k_max;
    Real        s1;

    A = sp_get(m,n,8);
    px = px_resize(px,n);
    MEM_STAT_REG(px,TYPE_PERM);
    for ( i = 0; i < A->m; i++ )
    {
        k_max = 1 + ((rand() >> 8) % 10);
        for ( k = 0; k < k_max; k++ )
        {
            j = (rand() >> 8) % A->n;
            s1 = rand()/((double)MAX_RAND);
            sp_set_val(A,i,j,s1);
        }
    }
    /* to make it likely that A is nonsingular, use pivot... */
    for ( i = 0; i < 2*A->n; i++ )
    {
        j = (rand() >> 8) % A->n;
        k = (rand() >> 8) % A->n;
        px_transp(px,j,k);
    }
    for ( i = 0; i < A->n; i++ )
        sp_set_val(A,i,px->pe[i],1.0);

    
    return A;
}
#endif

void stat_test1(par)
int par;
{
   static MAT *AT = MNULL;
   static VEC *xt1 = VNULL, *yt1 = VNULL;
   static VEC *xt2 = VNULL, *yt2 = VNULL;
   static VEC *xt3 = VNULL, *yt3 = VNULL;
   static VEC *xt4 = VNULL, *yt4 = VNULL;

   AT = m_resize(AT,10,10);
   xt1 = v_resize(xt1,10);
   yt1 = v_resize(yt1,10);
   xt2 = v_resize(xt2,10);
   yt2 = v_resize(yt2,10);
   xt3 = v_resize(xt3,10);
   yt3 = v_resize(yt3,10);
   xt4 = v_resize(xt4,10);
   yt4 = v_resize(yt4,10);

   MEM_STAT_REG(AT,TYPE_MAT);

#ifdef ANSI_C
   mem_stat_reg_vars(0,TYPE_VEC,__FILE__,__LINE__,&xt1,&xt2,&xt3,&xt4,&yt1,
		     &yt2,&yt3,&yt4,NULL);
#else
#ifdef VARARGS
   mem_stat_reg_vars(0,TYPE_VEC,__FILE__,__LINE__,&xt1,&xt2,&xt3,&xt4,&yt1,
		     &yt2,&yt3,&yt4,NULL);
#else
   MEM_STAT_REG(xt1,TYPE_VEC);
   MEM_STAT_REG(yt1,TYPE_VEC);
   MEM_STAT_REG(xt2,TYPE_VEC);
   MEM_STAT_REG(yt2,TYPE_VEC);
   MEM_STAT_REG(xt3,TYPE_VEC);
   MEM_STAT_REG(yt3,TYPE_VEC);
   MEM_STAT_REG(xt4,TYPE_VEC);
   MEM_STAT_REG(yt4,TYPE_VEC);
#endif
#endif

   v_rand(xt1);
   m_rand(AT);
   mv_mlt(AT,xt1,yt1);
   
}


void stat_test2(par)
int par;
{
   static PERM *px = PNULL;
   static IVEC *ixt = IVNULL, *iyt = IVNULL;
   
   px = px_resize(px,10);
   ixt = iv_resize(ixt,10);
   iyt = iv_resize(iyt,10);

   MEM_STAT_REG(px,TYPE_PERM);
   MEM_STAT_REG(ixt,TYPE_IVEC);
   MEM_STAT_REG(iyt,TYPE_IVEC);

   px_rand(px);
   px_inv(px,px);
}

#ifdef SPARSE
void stat_test3(par)
int par;
{
   static SPMAT *AT = (SPMAT *)NULL;
   static VEC *xt = VNULL, *yt = VNULL;
   static SPROW *r = (SPROW *) NULL;
   
   if (AT == (SPMAT *)NULL)
     AT = gen_non_symm(100,100);
   else
     AT = sp_resize(AT,100,100);
   xt = v_resize(xt,100);
   yt = v_resize(yt,100);
   if (r == NULL) r = sprow_get(100);

   MEM_STAT_REG(AT,TYPE_SPMAT);
   MEM_STAT_REG(xt,TYPE_VEC);
   MEM_STAT_REG(yt,TYPE_VEC);
   MEM_STAT_REG(r,TYPE_SPROW);

   v_rand(xt);
   sp_mv_mlt(AT,xt,yt);
   
}
#endif

#ifdef COMPLEX
void stat_test4(par)
int par;
{
   static ZMAT *AT = ZMNULL;
   static ZVEC *xt = ZVNULL, *yt = ZVNULL;
   
   AT = zm_resize(AT,10,10);
   xt = zv_resize(xt,10);
   yt = zv_resize(yt,10);

   MEM_STAT_REG(AT,TYPE_ZMAT);
   MEM_STAT_REG(xt,TYPE_ZVEC);
   MEM_STAT_REG(yt,TYPE_ZVEC);

   zv_rand(xt);
   zm_rand(AT);
   zmv_mlt(AT,xt,yt);
   
}
#endif


void main(argc, argv)
int     argc;
char    *argv[];
{
   VEC  *x = VNULL, *y = VNULL, *z = VNULL;
   PERM  *pi1 = PNULL, *pi2 = PNULL, *pi3 = PNULL;
   MAT   *A = MNULL, *B = MNULL, *C = MNULL;
#ifdef SPARSE
   SPMAT *sA, *sB;
   SPROW *r;
#endif
   IVEC *ix = IVNULL, *iy = IVNULL, *iz = IVNULL;
   int m,n,i,j,deg,k;
   Real s1,s2;
#ifdef COMPLEX
   ZVEC        *zx = ZVNULL, *zy = ZVNULL, *zz = ZVNULL;
   ZMAT        *zA = ZMNULL, *zB = ZMNULL, *zC = ZMNULL;
   complex     ONE;
#endif
   /* variables for testing attaching new lists of types  */
   FOO_1 *foo_1;
   FOO_2 *foo_2;


   mem_info_on(TRUE);

#if defined(ANSI_C) || defined(VARARGS)

   notice("vector initialize, copy & resize");
   
   n = v_get_vars(15,&x,&y,&z,(VEC **)NULL);
   if (n != 3) {
      errmesg("v_get_vars");
      printf(" n = %d (should be 3)\n",n);
   }

   v_rand(x);
   v_rand(y);