OS/2 Shareware BBS: 22 gnu/22-gnu.zip/mesch12a.zip/hsehldr.c

C/C++ Source or Header | 1994-01-13 | 5KB | 179 lines

/**************************************************************************
**
** 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.
**
***************************************************************************/


/*
        Files for matrix computations

    Householder transformation file. Contains routines for calculating
    householder transformations, applying them to vectors and matrices
    by both row & column.
*/

/* hsehldr.c 1.3 10/8/87 */
static    char    rcsid[] = "$Id: hsehldr.c,v 1.2 1994/01/13 05:36:29 des Exp $";

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


/* hhvec -- calulates Householder vector to eliminate all entries after the
    i0 entry of the vector vec. It is returned as out. May be in-situ */
VEC    *hhvec(vec,i0,beta,out,newval)
VEC    *vec,*out;
u_int    i0;
Real    *beta,*newval;
{
    Real    norm;

    out = _v_copy(vec,out,i0);
    norm = sqrt(_in_prod(out,out,i0));
    if ( norm <= 0.0 )
    {
        *beta = 0.0;
        return (out);
    }
    *beta = 1.0/(norm * (norm+fabs(out->ve[i0])));
    if ( out->ve[i0] > 0.0 )
        *newval = -norm;
    else
        *newval = norm;
    out->ve[i0] -= *newval;

    return (out);
}

/* hhtrvec -- apply Householder transformation to vector -- may be in-situ */
VEC    *hhtrvec(hh,beta,i0,in,out)
VEC    *hh,*in,*out;    /* hh = Householder vector */
u_int    i0;
double    beta;
{
    Real    scale;
    /* u_int    i; */

    if ( hh==(VEC *)NULL || in==(VEC *)NULL )
        error(E_NULL,"hhtrvec");
    if ( in->dim != hh->dim )
        error(E_SIZES,"hhtrvec");
    if ( i0 > in->dim )
        error(E_BOUNDS,"hhtrvec");

    scale = beta*_in_prod(hh,in,i0);
    out = v_copy(in,out);
    __mltadd__(&(out->ve[i0]),&(hh->ve[i0]),-scale,(int)(in->dim-i0));
    /************************************************************
    for ( i=i0; i<in->dim; i++ )
        out->ve[i] = in->ve[i] - scale*hh->ve[i];
    ************************************************************/

    return (out);
}

/* hhtrrows -- transform a matrix by a Householder vector by rows
    starting at row i0 from column j0 -- in-situ */
MAT    *hhtrrows(M,i0,j0,hh,beta)
MAT    *M;
u_int    i0, j0;
VEC    *hh;
double    beta;
{
    Real    ip, scale;
    int    i /*, j */;

    if ( M==(MAT *)NULL || hh==(VEC *)NULL )
        error(E_NULL,"hhtrrows");
    if ( M->n != hh->dim )
        error(E_RANGE,"hhtrrows");
    if ( i0 > M->m || j0 > M->n )
        error(E_BOUNDS,"hhtrrows");

    if ( beta == 0.0 )    return (M);

    /* for each row ... */
    for ( i = i0; i < M->m; i++ )
    {    /* compute inner product */
        ip = __ip__(&(M->me[i][j0]),&(hh->ve[j0]),(int)(M->n-j0));
        /**************************************************
        ip = 0.0;
        for ( j = j0; j < M->n; j++ )
            ip += M->me[i][j]*hh->ve[j];
        **************************************************/
        scale = beta*ip;
        if ( scale == 0.0 )
            continue;

        /* do operation */
        __mltadd__(&(M->me[i][j0]),&(hh->ve[j0]),-scale,
                            (int)(M->n-j0));
        /**************************************************
        for ( j = j0; j < M->n; j++ )
            M->me[i][j] -= scale*hh->ve[j];
        **************************************************/
    }

    return (M);
}


/* hhtrcols -- transform a matrix by a Householder vector by columns
    starting at row i0 from column j0 -- in-situ */
MAT    *hhtrcols(M,i0,j0,hh,beta)
MAT    *M;
u_int    i0, j0;
VEC    *hh;
double    beta;
{
    /* Real    ip, scale; */
    int    i /*, k */;
    static    VEC    *w = VNULL;

    if ( M==(MAT *)NULL || hh==(VEC *)NULL )
        error(E_NULL,"hhtrcols");
    if ( M->m != hh->dim )
        error(E_SIZES,"hhtrcols");
    if ( i0 > M->m || j0 > M->n )
        error(E_BOUNDS,"hhtrcols");

    if ( beta == 0.0 )    return (M);

    w = v_resize(w,M->n);
    MEM_STAT_REG(w,TYPE_VEC);
    v_zero(w);

    for ( i = i0; i < M->m; i++ )
        if ( hh->ve[i] != 0.0 )
        __mltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
                            (int)(M->n-j0));
    for ( i = i0; i < M->m; i++ )
        if ( hh->ve[i] != 0.0 )
        __mltadd__(&(M->me[i][j0]),&(w->ve[j0]),-beta*hh->ve[i],
                            (int)(M->n-j0));
    return (M);
}