ARM Club 1/ARM_CLUB_CD.iso/contents/apps/clib/progs/meschach/!Meschach/c/zmemory

Text File | 1994-04-04 | 14.9 KB | 645 lines

  va_end(ap);
   return i;
}



#elif VARARGS

#include <varargs.h>

/* To allocate memory to many arguments.
   The function should be called:
   v_get_vars(dim,&x,&y,&z,...,NULL);
   where
     int dim;
     ZVEC *x, *y, *z,...;
     The last argument should be NULL !
     dim is the length of vectors x,y,z,...
     returned value is equal to the number of allocated variables
     Other gec_... functions are similar.
*/

int zv_get_vars(va_alist) va_dcl
{
   va_list ap;
   int dim,i=0;
   ZVEC **par;

   va_start(ap);
   dim = va_arg(ap,int);
   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
      *par = zv_get(dim);
      i++;
   }

   va_end(ap);
   return i;
}



int zm_get_vars(va_alist) va_dcl
{
   va_list ap;
   int i=0, n, m;
   ZMAT **par;

   va_start(ap);
   m = va_arg(ap,int);
   n = va_arg(ap,int);
   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
      *par = zm_get(m,n);
      i++;
   }

   va_end(ap);
   return i;
}



/* To resize memory for many arguments.
   The function should be called:
   v_resize_vars(new_dim,&x,&y,&z,...,NULL);
   where
     int new_dim;
     ZVEC *x, *y, *z,...;
     The last argument should be NULL !
     rdim is the resized length of vectors x,y,z,...
     returned value is equal to the number of allocated variables.
     If one of x,y,z,.. arguments is NULL then memory is allocated to this
     argument.
     Other *_resize_list() functions are similar.
*/

int zv_resize_vars(va_alist) va_dcl
{
   va_list ap;
   int i=0, new_dim;
   ZVEC **par;

   va_start(ap);
   new_dim = va_arg(ap,int);
   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
      *par = zv_resize(*par,new_dim);
      i++;
   }

   va_end(ap);
   return i;
}


int zm_resize_vars(va_alist) va_dcl
{
   va_list ap;
   int i=0, m, n;
   ZMAT **par;

   va_start(ap);
   m = va_arg(ap,int);
   n = va_arg(ap,int);
   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
      *par = zm_resize(*par,m,n);
      i++;
   }

   va_end(ap);
   return i;
}



/* To deallocate memory for many arguments.
   The function should be called:
   v_free_vars(&x,&y,&z,...,NULL);
   where
     ZVEC *x, *y, *z,...;
     The last argument should be NULL !
     There must be at least one not NULL argument.
     returned value is equal to the number of allocated variables.
     Returned value of x,y,z,.. is VNULL.
     Other *_free_list() functions are similar.
*/

int zv_free_vars(va_alist) va_dcl
{
   va_list ap;
   int i=0;
   ZVEC **par;

   va_start(ap);
   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
      zv_free(*par);
      *par = ZVNULL;
      i++;
   }

   va_end(ap);
   return i;
}



int zm_free_vars(va_alist) va_dcl
{
   va_list ap;
   int i=0;
   ZMAT **par;

   va_start(ap);
   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
      zm_free(*par);
      *par = ZMNULL;
      i++;
   }

   va_end(ap);
   return i;
}


#endif

FileDataŵznormEÿÿÿ¬@9{
/**************************************************************************
**
** 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.
**
***************************************************************************/


/*
    A collection of functions for computing norms: scaled and unscaled
    Complex version
*/
static    char    rcsid[] = "$Id: znorm.c,v 1.1 1994/01/13 04:21:31 des Exp $";

#include    <stdio.h>
#include    <math.h>
#include    "zmatrix.h"



/* _zv_norm1 -- computes (scaled) 1-norms of vectors */
double    _zv_norm1(x,scale)
ZVEC    *x;
VEC    *scale;
{
    int    i, dim;
    Real    s, sum;

    if ( x == ZVNULL )
    error(E_NULL,"_zv_norm1");
    dim = x->dim;

    sum = 0.0;
    if ( scale == VNULL )
    for ( i = 0; i < dim; i++ )
        sum += zabs(x->ve[i]);
    else if ( scale->dim < dim )
    error(E_SIZES,"_zv_norm1");
    else
    for ( i = 0; i < dim; i++ )
    {
        s = scale->ve[i];
        sum += ( s== 0.0 ) ? zabs(x->ve[i]) : zabs(x->ve[i])/fabs(s);
    }

    return sum;
}

/* square -- returns x^2 */
/******************************
double    square(x)
double    x;
{    return x*x;    }
******************************/

#define    square(x)    ((x)*(x))

/* _zv_norm2 -- computes (scaled) 2-norm (Euclidean norm) of vectors */
double    _zv_norm2(x,scale)
ZVEC    *x;
VEC    *scale;
{
    int    i, dim;
    Real    s, sum;

    if ( x == ZVNULL )
    error(E_NULL,"_zv_norm2");
    dim = x->dim;

    sum = 0.0;
    if ( scale == VNULL )
    for ( i = 0; i < dim; i++ )
        sum += square(x->ve[i].re) + square(x->ve[i].im);
    else if ( scale->dim < dim )
    error(E_SIZES,"_v_norm2");
    else
    for ( i = 0; i < dim; i++ )
    {
        s = scale->ve[i];
        sum += ( s== 0.0 ) ? square(x->ve[i].re) + square(x->ve[i].im) :
        (square(x->ve[i].re) + square(x->ve[i].im))/square(s);
    }

    return sqrt(sum);
}

#define    max(a,b)    ((a) > (b) ? (a) : (b))

/* _zv_norm_inf -- computes (scaled) infinity-norm (supremum norm) of vectors */
double    _zv_norm_inf(x,scale)
ZVEC    *x;
VEC    *scale;
{
    int    i, dim;
    Real    s, maxval, tmp;

    if ( x == ZVNULL )
    error(E_NULL,"_zv_norm_inf");
    dim = x->dim;

    maxval = 0.0;
    if ( scale == VNULL )
    for ( i = 0; i < dim; i++ )
    {
        tmp = zabs(x->ve[i]);
        maxval = max(maxval,tmp);
    }
    else if ( scale->dim < dim )
    error(E_SIZES,"_zv_norm_inf");
    else
    for ( i = 0; i < dim; i++ )
    {
        s = scale->ve[i];
        tmp = ( s == 0.0 ) ? zabs(x->ve[i]) : zabs(x->ve[i])/fabs(s);
        maxval = max(maxval,tmp);
    }

    return maxval;
}

/* zm_norm1 -- compute matrix 1-norm -- unscaled
    -- complex version */
double    zm_norm1(A)
ZMAT    *A;
{
    int    i, j, m, n;
    Real    maxval, sum;

    if ( A == ZMNULL )
    error(E_NULL,"zm_norm1");

    m = A->m;    n = A->n;
    maxval = 0.0;

    for ( j = 0; j < n; j++ )
    {
    sum = 0.0;
    for ( i = 0; i < m; i ++ )
        sum += zabs(A->me[i][j]);
    maxval = max(maxval,sum);
    }

    return maxval;
}

/* zm_norm_inf -- compute matrix infinity-norm -- unscaled
    -- complex version */
double    zm_norm_inf(A)
ZMAT    *A;
{
    int    i, j, m, n;
    Real    maxval, sum;

    if ( A == ZMNULL )
    error(E_NULL,"zm_norm_inf");

    m = A->m;    n = A->n;
    maxval = 0.0;

    for ( i = 0; i < m; i++ )
    {
    sum = 0.0;
    for ( j = 0; j < n; j ++ )
        sum += zabs(A->me[i][j]);
    maxval = max(maxval,sum);
    }

    return maxval;
}

/* zm_norm_frob -- compute matrix frobenius-norm -- unscaled */
double    zm_norm_frob(A)
ZMAT    *A;
{
    int    i, j, m, n;
    Real    sum;

    if ( A == ZMNULL )
    error(E_NULL,"zm_norm_frob");

    m = A->m;    n = A->n;
    sum = 0.0;

    for ( i = 0; i < m; i++ )
    for ( j = 0; j < n; j ++ )
        sum += square(A->me[i][j].re) + square(A->me[i][j].im);

    return sqrt(sum);
}

FileDataŵzqrfctrÜ5Eÿÿÿ¬Þøÿ=
/**************************************************************************
**
** 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.
**
***************************************************************************/

/*
  This file contains the routines needed to perform QR factorisation
  of matrices, as well as Householder transformations.
  The internal "factored form" of a matrix A is not quite standard.
  The diagonal of A is replaced by the diagonal of R -- not by the 1st non-zero
  entries of the Householder vectors. The 1st non-zero entries are held in
  the diag parameter of QRfactor(). The reason for this non-standard
  representation is that it enables direct use of the Usolve() function
  rather than requiring that  a seperate function be written just for this case.
  See, e.g., QRsolve() below for more details.

  Complex version

*/

static    char    rcsid[] = "$Id: zqrfctr.c,v 1.1 1994/01/13 04:21:22 des Exp $";

#include    <stdio.h>
#include    <math.h>
#include    "zmatrix.h"
#include    "zmatrix2.h"


#define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)


#define        sign(x)    ((x) > 0.0 ? 1 : ((x) < 0.0 ? -1 : 0 ))

/* Note: The usual representation of a Householder transformation is taken
   to be:
   P = I - beta.u.u*
   where beta = 2/(u*.u) and u is called the Householder vector
   (u* is the conjugate transposed vector of u
*/

/* zQRfactor -- forms the QR factorisation of A
    -- factorisation stored in compact form as described above
    (not quite standard format) */
ZMAT    *zQRfactor(A,diag)
ZMAT    *A;
ZVEC    *diag;
{
    u_int    k,limit;
    Real    beta;
    static    ZVEC    *tmp1=ZVNULL;

    if ( ! A || ! diag )
    error(E_NULL,"zQRfactor");
    limit = min(A->m,A->n);
    if ( diag->dim < limit )
    error(E_SIZES,"zQRfactor");

    tmp1 = zv_resize(tmp1,A->m);
    MEM_STAT_REG(tmp1,TYPE_ZVEC);

    for ( k=0; k<limit; k++ )
    {
    /* get H/holder vector for the k-th column */
    zget_col(A,k,tmp1);
    /* hhvec(tmp1,k,&beta->ve[k],tmp1,&A->me[k][k]); */
    zhhvec(tmp1,k,&beta,tmp1,&A->me[k][k]);
    diag->ve[k] = tmp1->ve[k];

    /* apply H/holder vector to remaining columns */
    /* hhtrcols(A,k,k+1,tmp1,beta->ve[k]); */
    tracecatch(zhhtrcols(A,k,k+1,tmp1,beta),"zQRfactor");
    }

    return (A);
}

/* zQRCPfactor -- forms the QR factorisation of A with column pivoting
   -- factorisation stored in compact form as described above
   ( not quite standard format )                */
ZMAT    *zQRCPfactor(A,diag,px)
ZMAT    *A;
ZVEC    *diag;
PERM    *px;
{
    u_int    i, i_max, j, k, limit;
    static    ZVEC    *tmp1=ZVNULL, *tmp2=ZVNULL;
    static    VEC    *gamma=VNULL;
    Real     beta;
    Real    maxgamma, sum, tmp;
    complex    ztmp;

    if ( ! A || ! diag || ! px )
    error(E_NULL,"QRCPfactor");
    limit = min(A->m,A->n);
    if ( diag->dim < limit || px->size != A->n )
    error(E_SIZES,"QRCPfactor");

    tmp1 = zv_resize(tmp1,A->m);
    tmp2 = zv_resize(tmp2,A->m);
    gamma = v_resize(gamma,A->n);
    MEM_STAT_REG(tmp1,TYPE_ZVEC);
    MEM_STAT_REG(tmp2,TYPE_ZVEC);
    MEM_STAT_REG(gamma,TYPE_VEC);

    /* initialise gamma and px */
    for ( j=0; j<A->n; j++ )
    {
    px->pe[j] = j;
    sum = 0.0;
    for ( i=0; i<A->m; i++ )
        sum += square(A->me[i][j].re) + square(A->me[i][j].im);
    gamma->ve[j] = sum;
    }

    for ( k=0; k<limit; k++ )
    {
    /* find "best" column to use */
    i_max = k;    maxgamma = gamma->ve[k];
    for ( i=k+1; i<A->n; i++ )
        /* Loop invariant:maxgamma=gamma[i_max]
           >=gamma[l];l=k,...,i-1 */
        if ( gamma->ve[i] > maxgamma )
        {    maxgamma = gamma->ve[i]; i_max = i;    }

    /* swap columns if necessary */
    if ( i_max != k )
    {
        /* swap gamma values */
        tmp = gamma->ve[k];
        gamma->ve[k] = gamma->ve[i_max];
        gamma->ve[i_max] = tmp;

        /* update column permutation */
        px_transp(px,k,i_max);

        /* swap columns of A */
        for ( i=0; i<A->m; i++ )
        {
        ztmp = A->me[i][k];
        A->me[i][k] = A->me[i][i_max];
        A->me[i][i_max] = ztmp;
        }
    }

    /* get H/holder vector for the k-th column */
    zget_col(A,k,tmp1);
    /* hhvec(tmp1,k,&beta->ve[k],tmp1,&A->me[k][k]); */
    zhhvec(tmp1,k,&beta,tmp1,&A->me[k][k]);
    diag->ve[k] = tmp1->ve[k];

    /* apply H/holder vector to remaining columns */
    /* hhtrcols(A,k,k+1,tmp1,beta->ve[k]); */
    zhhtrcols(A,k,k+1,tmp1,beta);

    /* update gamma values */
    for ( j=k+1; j<A->n; j++ )
        gamma->ve[j] -= square(A->me[k][j].re)+square(A->me[k][j].im);
    }

    return (A);
}

/* zQsolve -- solves Qx = b, Q is an orthogonal matrix stored in compact
    form a la QRfactor()
    -- may be in-situ */
ZVEC    *_zQsolve(QR,diag,b,x,tmp)
ZMAT    *QR;
ZVEC    *diag, *b, *x, *tmp;
{
    u_int    dynamic;
    int        k, limit;
    Real    beta, r_ii, tmp_val;

    limit = min(QR->m,QR->n);
    dynamic = FALSE;
    if ( ! QR || ! diag || ! b )
    error(E_NULL,"_zQsolve");
    if ( diag->dim < limit || b->dim != QR->m )
    error(E_SIZES,"_zQsolve");
    x = zv_resize(x,QR->m);
    if ( tmp == ZVNULL )
    dynamic = TRUE;
    tmp = zv_resize(tmp,QR->m);

    /* apply H/holder transforms in normal order */
    x = zv_copy(b,x);
    for ( k = 0 ; k < limit ; k++ )
    {
    zget_col(QR,k,tmp);
    r_ii = zabs(tmp->ve[k]);
    tmp->ve[k] = diag->ve[k];
    tmp_val = (r_ii*zabs(diag->ve[k]));
    beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
    /* hhtrvec(tmp,beta->ve[k],k,x,x); */
    zhhtrvec(tmp,beta,k,x,x);
    }

    if ( dynamic )
    ZV_FREE(tmp);

    return (x);
}

/* zmakeQ -- constructs orthogonal matrix from Householder vectors stored in
   compact QR form */
ZMAT    *zmakeQ(QR,diag,Qout)
ZMAT    *QR,*Qout;
ZVEC    *diag;
{
    static    ZVEC    *tmp1=ZVNULL,*tmp2=ZVNULL;
    u_int    i, limit;
    Real    beta, r_ii, tmp_val;
    int    j;

    limit = min(QR->m,QR->n);
    if ( ! QR || ! diag )
    error(E_NULL,"zmakeQ");
    if ( diag->dim < limit )
    error(E_SIZES,"zmakeQ");
    Qout = zm_resize(Qout,QR->m,QR->m);

    tmp1 = zv_resize(tmp1,QR->m);    /* contains basis vec & columns of Q */
    tmp2 = zv_resize(tmp2,QR->m);    /* contains H/holder vectors */
    MEM_STAT_REG(tmp1,TYPE_ZVEC);
    MEM_STAT_REG(tmp2,TYPE_ZVEC);

    for ( i=0; i<QR->m ; i++ )
    {    /* get i-th column of Q */
    /* set up tmp1 as i-th basis vector */
    for ( j=0; j<QR->m ; j++ )
        tmp1->ve[j].re = tmp1->ve[j].im = 0.0;
    tmp1->ve[i].re = 1.0;

    /* apply H/h transforms in reverse order */
    for ( j=limit-1; j>=0; j-- )
    {
        zget_col(QR,j,tmp2);
        r_ii = zabs(tmp2->ve[j]);
        tmp2->ve[j] = diag->ve[j];
        tmp_val = (r_ii*zabs(diag->ve[j]));
        beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
        /* hhtrvec(tmp2,beta->ve[j],j,tmp1,tmp1); */
        zhhtrvec(tmp2,beta,j,tmp1,tmp1);
    }

    /* insert into Q */
    zset_col(Qout,i,tmp1);
    }

    return (Qout);
}

/* zmakeR -- constructs upper triangu