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

Text File | 1994-03-08 | 23.2 KB | 1,024 lines

05:34:35 des Exp $";

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


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

    if ( x == (VEC *)NULL )
        error(E_NULL,"_v_norm1");
    dim = x->dim;

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

    return sum;
}

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

/* cube -- returns x^3 */
double cube(x)
double x;
{  return x*x*x;   }

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

    if ( x == (VEC *)NULL )
        error(E_NULL,"_v_norm2");
    dim = x->dim;

    sum = 0.0;
    if ( scale == (VEC *)NULL )
        for ( i = 0; i < dim; i++ )
            sum += square(x->ve[i]);
    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]) :
                            square(x->ve[i]/s);
        }

    return sqrt(sum);
}

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

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

    if ( x == (VEC *)NULL )
        error(E_NULL,"_v_norm_inf");
    dim = x->dim;

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

    return maxval;
}

/* m_norm1 -- compute matrix 1-norm -- unscaled */
double    m_norm1(A)
MAT    *A;
{
    int    i, j, m, n;
    Real    maxval, sum;

    if ( A == (MAT *)NULL )
        error(E_NULL,"m_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 += fabs(A->me[i][j]);
        maxval = max(maxval,sum);
    }

    return maxval;
}

/* m_norm_inf -- compute matrix infinity-norm -- unscaled */
double    m_norm_inf(A)
MAT    *A;
{
    int    i, j, m, n;
    Real    maxval, sum;

    if ( A == (MAT *)NULL )
        error(E_NULL,"m_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 += fabs(A->me[i][j]);
        maxval = max(maxval,sum);
    }

    return maxval;
}

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

    if ( A == (MAT *)NULL )
        error(E_NULL,"m_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]);

    return sqrt(sum);
}

FileDataŵotherioŵEÿÿÿ‘ñ?ñã
/**************************************************************************
**
** 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.
**
***************************************************************************/


/*
    File for doing assorted I/O operations not invlolving
    MAT/VEC/PERM objects
*/
static    char    rcsid[] = "$Id: otherio.c,v 1.2 1994/01/13 05:34:52 des Exp $";

#include    <stdio.h>
#include    <ctype.h>
#include    "matrix.h"



/* scratch area -- enough for a single line */
static    char    scratch[MAXLINE+1];

/* default value for fy_or_n */
static    int    y_n_dflt = TRUE;

/* fy_or_n -- yes-or-no to question is string s
    -- question written to stderr, input from fp
    -- if fp is NOT a tty then return y_n_dflt */
int    fy_or_n(fp,s)
FILE    *fp;
char    *s;
{
    char    *cp;

    if ( ! isatty(fileno(fp)) )
        return y_n_dflt;

    for ( ; ; )
    {
        fprintf(stderr,"%s (y/n) ? ",s);
        if ( fgets(scratch,MAXLINE,fp)==NULL )
            error(E_INPUT,"fy_or_n");
        cp = scratch;
        while ( isspace(*cp) )
            cp++;
        if ( *cp == 'y' || *cp == 'Y' )
            return TRUE;
        if ( *cp == 'n' || *cp == 'N' )
            return FALSE;
        fprintf(stderr,"Please reply with 'y' or 'Y' for yes ");
        fprintf(stderr,"and 'n' or 'N' for no.\n");
    }
}

/* yn_dflt -- sets the value of y_n_dflt to val */
int    yn_dflt(val)
int    val;
{    return y_n_dflt = val;        }

/* fin_int -- return integer read from file/stream fp
    -- prompt s on stderr if fp is a tty
    -- check that x lies between low and high: re-prompt if
        fp is a tty, error exit otherwise
    -- ignore check if low > high        */
int    fin_int(fp,s,low,high)
FILE    *fp;
char    *s;
int    low, high;
{
    int    retcode, x;

    if ( ! isatty(fileno(fp)) )
    {
        skipjunk(fp);
        if ( (retcode=fscanf(fp,"%d",&x)) == EOF )
            error(E_INPUT,"fin_int");
        if ( retcode <= 0 )
            error(E_FORMAT,"fin_int");
        if ( low <= high && ( x < low || x > high ) )
            error(E_BOUNDS,"fin_int");
        return x;
    }

    for ( ; ; )
    {
        fprintf(stderr,"%s: ",s);
        if ( fgets(scratch,MAXLINE,stdin)==NULL )
            error(E_INPUT,"fin_int");
        retcode = sscanf(scratch,"%d",&x);
        if ( ( retcode==1 && low > high ) ||
                    ( x >= low && x <= high ) )
            return x;
        fprintf(stderr,"Please type an integer in range [%d,%d].\n",
                            low,high);
    }
}


/* fin_double -- return double read from file/stream fp
    -- prompt s on stderr if fp is a tty
    -- check that x lies between low and high: re-prompt if
        fp is a tty, error exit otherwise
    -- ignore check if low > high        */
double    fin_double(fp,s,low,high)
FILE    *fp;
char    *s;
double    low, high;
{
    Real    retcode, x;

    if ( ! isatty(fileno(fp)) )
    {
        skipjunk(fp);
#if REAL == DOUBLE
        if ( (retcode=fscanf(fp,"%lf",&x)) == EOF )
#elif REAL == FLOAT
        if ( (retcode=fscanf(fp,"%f",&x)) == EOF )
#endif
            error(E_INPUT,"fin_double");
        if ( retcode <= 0 )
            error(E_FORMAT,"fin_double");
        if ( low <= high && ( x < low || x > high ) )
            error(E_BOUNDS,"fin_double");
        return (double)x;
    }

    for ( ; ; )
    {
        fprintf(stderr,"%s: ",s);
        if ( fgets(scratch,MAXLINE,stdin)==NULL )
            error(E_INPUT,"fin_double");
#if REAL == DOUBLE
        retcode = sscanf(scratch,"%lf",&x);
#elif REAL == FLOAT
        retcode = sscanf(scratch,"%f",&x);
#endif
        if ( ( retcode==1 && low > high ) ||
                    ( x >= low && x <= high ) )
            return (double)x;
        fprintf(stderr,"Please type an double in range [%g,%g].\n",
                            low,high);
    }
}


FileDataŵpxopIEÿÿÿàÆ-;èº
/**************************************************************************
**
** 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.
**
***************************************************************************/


/* pxop.c 1.5 12/03/87 */


#include    <stdio.h>
#include    "matrix.h"

static    char    rcsid[] = "$Id: pxop.c,v 1.5 1994/03/23 23:58:50 des Exp $";

/**********************************************************************
Note: A permutation is often interpreted as a matrix
        (i.e. a permutation matrix).
    A permutation px represents a permutation matrix P where
        P[i][j] == 1 if and only if px->pe[i] == j
**********************************************************************/


/* px_inv -- invert permutation -- in situ
    -- taken from ACM Collected Algorithms #250 */
PERM    *px_inv(px,out)
PERM    *px, *out;
{
    int    i, j, k, n, *p;

    out = px_copy(px, out);
    n = out->size;
    p = (int *)(out->pe);
    for ( n--; n>=0; n-- )
    {
    i = p[n];
    if ( i < 0 )    p[n] = -1 - i;
    else if ( i != n )
    {
        k = n;
        while (TRUE)
        {
        if ( i < 0 || i >= out->size )
            error(E_BOUNDS,"px_inv");
        j = p[i];    p[i] = -1 - k;
        if ( j == n )
        {    p[n] = i;    break;        }
        k = i;        i = j;
        }
    }
    }
    return out;
}

/* px_mlt -- permutation multiplication (composition) */
PERM    *px_mlt(px1,px2,out)
PERM    *px1,*px2,*out;
{
    u_int    i,size;

    if ( px1==(PERM *)NULL || px2==(PERM *)NULL )
    error(E_NULL,"px_mlt");
    if ( px1->size != px2->size )
    error(E_SIZES,"px_mlt");
    if ( px1 == out || px2 == out )
    error(E_INSITU,"px_mlt");
    if ( out==(PERM *)NULL || out->size < px1->size )
    out = px_resize(out,px1->size);

    size = px1->size;
    for ( i=0; i<size; i++ )
    if ( px2->pe[i] >= size )
        error(E_BOUNDS,"px_mlt");
    else
        out->pe[i] = px1->pe[px2->pe[i]];

    return out;
}

/* px_vec -- permute vector */
VEC    *px_vec(px,vector,out)
PERM    *px;
VEC    *vector,*out;
{
    u_int    old_i, i, size, start;
    Real    tmp;

    if ( px==(PERM *)NULL || vector==(VEC *)NULL )
    error(E_NULL,"px_vec");
    if ( px->size > vector->dim )
    error(E_SIZES,"px_vec");
    if ( out==(VEC *)NULL || out->dim < vector->dim )
    out = v_resize(out,vector->dim);

    size = px->size;
    if ( size == 0 )
    return v_copy(vector,out);
    if ( out != vector )
    {
    for ( i=0; i<size; i++ )
        if ( px->pe[i] >= size )
        error(E_BOUNDS,"px_vec");
        else
        out->ve[i] = vector->ve[px->pe[i]];
    }
    else
    {    /* in situ algorithm */
    start = 0;
    while ( start < size )
    {
        old_i = start;
        i = px->pe[old_i];
        if ( i >= size )
        {
        start++;
        continue;
        }
        tmp = vector->ve[start];
        while ( TRUE )
        {
        vector->ve[old_i] = vector->ve[i];
        px->pe[old_i] = i+size;
        old_i = i;
        i = px->pe[old_i];
        if ( i >= size )
            break;
        if ( i == start )
        {
            vector->ve[old_i] = tmp;
            px->pe[old_i] = i+size;
            break;
        }
        }
        start++;
    }

    for ( i = 0; i < size; i++ )
        if ( px->pe[i] < size )
        error(E_BOUNDS,"px_vec");
        else
        px->pe[i] = px->pe[i]-size;
    }

    return out;
}

/* pxinv_vec -- apply the inverse of px to x, returning the result in out */
VEC    *pxinv_vec(px,x,out)
PERM    *px;
VEC    *x, *out;
{
    u_int    i, size;

    if ( ! px || ! x )
    error(E_NULL,"pxinv_vec");
    if ( px->size > x->dim )
    error(E_SIZES,"pxinv_vec");
    /* if ( x == out )
    error(E_INSITU,"pxinv_vec"); */
    if ( ! out || out->dim < x->dim )
    out = v_resize(out,x->dim);

    size = px->size;
    if ( size == 0 )
    return v_copy(x,out);
    if ( out != x )
    {
    for ( i=0; i<size; i++ )
        if ( px->pe[i] >= size )
        error(E_BOUNDS,"pxinv_vec");
        else
        out->ve[px->pe[i]] = x->ve[i];
    }
    else
    {    /* in situ algorithm --- cheat's way out */
    px_inv(px,px);
    px_vec(px,x,out);
    px_inv(px,px);
    }

    return out;
}



/* px_transp -- transpose elements of permutation
        -- Really multiplying a permutation by a transposition */
PERM    *px_transp(px,i1,i2)
PERM    *px;        /* permutation to transpose */
u_int    i1,i2;        /* elements to transpose */
{
    u_int    temp;

    if ( px==(PERM *)NULL )
        error(E_NULL,"px_transp");

    if ( i1 < px->size && i2 < px->size )
    {
        temp = px->pe[i1];
        px->pe[i1] = px->pe[i2];
        px->pe[i2] = temp;
    }

    return px;
}

/* myqsort -- a cheap implementation of Quicksort on integers
        -- returns number of swaps */
static int myqsort(a,num)
int    *a, num;
{
    int    i, j, tmp, v;
    int    numswaps;

    numswaps = 0;
    if ( num <= 1 )
        return 0;

    i = 0;    j = num;    v = a[0];
    for ( ; ; )
    {
        while ( a[++i] < v )
            ;
        while ( a[--j] > v )
            ;
        if ( i >= j )    break;

        tmp = a[i];
        a[i] = a[j];
        a[j] = tmp;
        numswaps++;
    }

    tmp = a[0];
    a[0] = a[j];
    a[j] = tmp;
    if ( j != 0 )
        numswaps++;

    numswaps += myqsort(&a[0],j);
    numswaps += myqsort(&a[j+1],num-(j+1));

    return numswaps;
}


/* px_sign -- compute the ``sign'' of a permutation = +/-1 where
        px is the product of an even/odd # transpositions */
int    px_sign(px)
PERM    *px;
{
    int    numtransp;
    PERM    *px2;

    if ( px==(PERM *)NULL )
        error(E_NULL,"px_sign");
    px2 = px_copy(px,PNULL);
    numtransp = myqsort(px2->pe,px2->size);
    px_free(px2);

    return ( numtransp % 2 ) ? -1 : 1;
}


/* px_cols -- permute columns of matrix A; out = A.px'
    -- May NOT be in situ */
MAT    *px_cols(px,A,out)
PERM    *px;
MAT    *A, *out;
{
    int    i, j, m, n, px_j;
    Real    **A_me, **out_me;
#ifdef ANSI_C
    MAT    *m_get(int, int);
#else
    extern MAT    *m_get();
#endif

    if ( ! A || ! px )
        error(E_NULL,"px_cols");
    if ( px->size != A->n )
        error(E_SIZES,"px_cols");
    if ( A == out )
        error(E_INSITU,"px_cols");
    m = A->m;    n = A->n;
    if ( ! out || out->m != m || out->n != n )
        out = m_get(m,n);
    A_me = A->me;    out_me = out->me;

    for ( j = 0; j < n; j++ )
    {
        px_j = px->pe[j];
        if ( px_j >= n )
            error(E_BOUNDS,"px_cols");
        for ( i = 0; i < m; i++ )
            out_me[i][px_j] = A_me[i][j];
    }

    return out;
}

/* px_rows -- permute columns of matrix A; out = px.A
    -- May NOT be in situ */
MAT    *px_rows(px,A,out)
PERM    *px;
MAT    *A, *out;
{
    int    i, j, m, n, px_i;
    Real    **A_me, **out_me;
#ifdef ANSI_C
    MAT    *m_get(int, int);
#else
    extern MAT    *m_get();
#endif

    if ( ! A || ! px )
        error(E_NULL,"px_rows");
    if ( px->size != A->m )
        error(E_SIZES,"px_rows");
    if ( A == out )
        error(E_INSITU,"px_rows");
    m = A->m;    n = A->n;
    if ( ! out || out->m != m || out->n != n )
        out = m_get(m,n);
    A_me = A->me;    out_me = out->me;

    for ( i = 0; i < m; i++ )
    {
        px_i = px->pe[i];
        if ( px_i >= m )
            error(E_BOUNDS,"px_rows");
        for ( j = 0; j < n; j++ )
            out_me[i][j] = A_me[px_i][j];
    }

    return out;
}

FileDataŵqrfactorE4EÿÿÿØþ?´!
/**************************************************************************
**
** 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.

*/


static    char    rcsid[] = "$Id: qrfactor.c,v 1.5 1994/01/13 05:35:07 des Exp $";

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





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

extern    VEC    *Usolve();    /* See matrix2.h */

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

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

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

    tmp1 = v_resize(tmp1,A->m);
    MEM_STAT_REG(tmp1,TYPE_VEC);

    for ( k=0; k<limit; k++ )
    {
    /* get H/holder vector for the k-th column */
    get_col(A,k,tmp1);
    /* hhvec(tmp1,k,&beta->ve[k],tmp1,&A->me[k][k]); */
    hhvec(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]); */
    hhtrcols(A,k,k+1,tmp1,beta);
    }

    return (A);
}

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

    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 = v_resize(tmp1,A->m);
    tmp2 = v_resize(tmp2,A->m);
    gamma = v_resize(gamma,A->n);
    MEM_STAT_REG(tmp1,TYPE_VEC);
    MEM_STAT_REG(tmp2,TYPE_VEC);
    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]);
    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++ )
        {
        tmp = A->me[i][k];
        A->me[i][k] = A->me[i][i_max];
        A->me[i][i_max] = tmp;
        }
    }

    /* get H/holder vector for the k-th column */
    get_col(A,k,tmp1);
    /* hhvec(tmp1,k,&beta->ve[k],tmp1,&A->me[k][k]); */
    hhvec(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]); */
    hhtrcols(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]);
    }

    return (A);
}

/* Qsolve -- solves Qx = b, Q is an orthogonal matrix stored in compact
   form a la QRfactor() -- may be in-situ */
/* VEC    *_Qsolve(QR,diag,beta,b,x,tmp) */
VEC    *_Qsolve(QR,diag,b,x,tmp)
MAT    *QR;
VEC    *diag /* ,*beta */ , *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,"_Qsolve");
    if ( diag->dim < limit || b->dim != QR->m )
    error(E_SIZES,"_Qsolve");
    x = v_resize(x,QR->m);
    if ( tmp == VNULL )
    dynamic = TRUE;
    tmption 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.
**
***************************************************************************/


/*
  Sparse matrix Bunch--Kaufman--Parlett factorisation and solve
  Radical revision started Thu 05th Nov 1992, 09:36:12 AM
  to use Karen George's suggestion of leaving the the row elements unordered
  Radical revision completed Mon 07th Dec 1992, 10:59:57 AM
*/

static    char    rcsid[] = "$Id: spbkp.c,v 1.5 1994/01/13 05:44:35 des Exp $";

#include    <stdio.h>
#include    <math.h>
#include    "matrix.h"
#include    "sparse.h"
#include        "sparse2.h"


#ifdef MALLOCDECL
#include <malloc.h>
#endif

#define alpha    0.6403882032022076 /* = (1+sqrt(17))/8 */


#define    btos(x)    ((x) ? "TRUE" : "FALSE")

/* assume no use of sqr() uses side-effects */
#define    sqr(x)    ((x)*(x))

/* unord_get_idx -- returns index (encoded if entry not allocated)
    of the element of row r with column j
    -- uses linear search */
int    unord_get_idx(r,j)
SPROW    *r;
int    j;
{
    int        idx;
    row_elt    *e;

    if ( ! r || ! r->elt )
    error(E_NULL,"unord_get_idx");
    for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
    if ( e->col == j )
        break;
    if ( idx >= r->len )
    return -(r->len+2);
    else
    return idx;
}

/* unord_get_val -- returns value of the (i,j) entry of A
    -- same assumptions as unord_get_idx() */
double    unord_get_val(A,i,j)
SPMAT    *A;
int    i, j;
{
    SPROW    *r;
    int        idx;

    if ( ! A )
    error(E_NULL,"unord_get_val");
    if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
    error(E_BOUNDS,"unord_get_val");

    r = &(A->row[i]);
    idx = unord_get_idx(r,j);
    if ( idx < 0 )
    return 0.0;
    else
    return r->elt[idx].val;
}


/* bkp_swap_elt -- swaps the (i,j) with the (k,l) entry of sparse matrix
    -- either or both of the entries may be unallocated */
static SPMAT    *bkp_swap_elt(A,i1,j1,idx1,i2,j2,idx2)
SPMAT    *A;
int    i1, j1, idx1, i2, j2, idx2;
{
    int        tmp_row, tmp_idx;
    SPROW    *r1, *r2;
    row_elt    *e1, *e2;
    Real    tmp;

    if ( ! A )
    error(E_NULL,"bkp_swap_elt");

    if ( i1 < 0 || j1 < 0 || i2 < 0 || j2 < 0 ||
     i1 >= A->m || j1 >= A->n || i2 >= A->m || j2 >= A->n )
    {
    error(E_BOUNDS,"bkp_swap_elt");
    }

    if ( i1 == i2 && j1 == j2 )
    return A;
    if ( idx1 < 0 && idx2 < 0 )    /* neither allocated */
    return A;

    r1 = &(A->row[i1]);        r2 = &(A->row[i2]);
    /* if ( idx1 >= r1->len || idx2 >= r2->len )
    error(E_BOUNDS,"bkp_swap_elt"); */
    if ( idx1 < 0 )    /* assume not allocated */
    {
    idx1 = r1->len;
    if ( idx1 >= r1->maxlen )
    {    tracecatch(sprow_xpd(r1,2*r1->maxlen+1,TYPE_SPMAT),
            "bkp_swap_elt");    }
    r1->len = idx1+1;
    r1->elt[idx1].col = j1;
    r1->elt[idx1].val = 0.0;
    /* now patch up column access path */
    tmp_row = -1;    tmp_idx = j1;
    chase_col(A,j1,&tmp_row,&tmp_idx,i1-1);

    if ( tmp_row < 0 )
    {
        r1->elt[idx1].nxt_row = A->start_row[j1];
        r1->elt[idx1].nxt_idx = A->start_idx[j1];
        A->start_row[j1] = i1;
        A->start_idx[j1] = idx1;
    }
    else
    {
        row_elt    *tmp_e;

        tmp_e = &(A->