zhsehldr.c

Go to the documentation of this file.

#include <../../nrnconf.h>
 
/**************************************************************************
**
** 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.
 
    Complex version
*/
 
static  char    rcsid[] = "zhsehldr.c,v 1.1 1997/12/04 17:56:08 hines Exp";
 
#include    <stdio.h>
#include    "zmatrix.h"
#include        "zmatrix2.h"
#include    <math.h>
 
#define is_zero(z)  ((z).re == 0.0 && (z).im == 0.0)
 
/* zhhvec -- calulates Householder vector to eliminate all entries after the
    i0 entry of the vector vec. It is returned as out. May be in-situ */
ZVEC    *zhhvec(vec,i0,beta,out,newval)
ZVEC    *vec,*out;
int i0;
Real    *beta;
complex *newval;
{
    complex tmp;
    Real    norm, abs_val;
 
    if ( i0 < 0 || i0 >= vec->dim )
        error(E_BOUNDS,"zhhvec");
    out = _zv_copy(vec,out,i0);
    tmp = _zin_prod(out,out,i0,Z_CONJ);
    if ( tmp.re <= 0.0 )
    {
        *beta = 0.0;
        *newval = out->ve[i0];
        return (out);
    }
    norm = sqrt(tmp.re);
    abs_val = zabs(out->ve[i0]);
    *beta = 1.0/(norm * (norm+abs_val));
    if ( abs_val == 0.0 )
    {
      newval->re = norm;
      newval->im = 0.0;
    }
    else
    { 
      abs_val = -norm / abs_val;
      newval->re = abs_val*out->ve[i0].re;
      newval->im = abs_val*out->ve[i0].im;
    }   abs_val = -norm / abs_val;
    out->ve[i0].re -= newval->re;
    out->ve[i0].im -= newval->im;
 
    return (out);
}
 
/* zhhtrvec -- apply Householder transformation to vector -- may be in-situ */
ZVEC    *zhhtrvec(hh,beta,i0,in,out)
ZVEC    *hh,*in,*out;   /* hh = Householder vector */
int i0;
double  beta;
{
    complex scale, tmp;
    /* u_int    i; */
 
    if ( hh==ZVNULL || in==ZVNULL )
        error(E_NULL,"zhhtrvec");
    if ( in->dim != hh->dim )
        error(E_SIZES,"zhhtrvec");
    if ( i0 < 0 || i0 > in->dim )
        error(E_BOUNDS,"zhhvec");
 
    tmp = _zin_prod(hh,in,i0,Z_CONJ);
    scale.re = -beta*tmp.re;
    scale.im = -beta*tmp.im;
    out = zv_copy(in,out);
    __zmltadd__(&(out->ve[i0]),&(hh->ve[i0]),scale,
            (int)(in->dim-i0),Z_NOCONJ);
    /************************************************************
    for ( i=i0; i<in->dim; i++ )
        out->ve[i] = in->ve[i] - scale*hh->ve[i];
    ************************************************************/
 
    return (out);
}
 
/* zhhtrrows -- transform a matrix by a Householder vector by rows
    starting at row i0 from column j0 -- in-situ */
ZMAT    *zhhtrrows(M,i0,j0,hh,beta)
ZMAT    *M;
int i0, j0;
ZVEC    *hh;
double  beta;
{
    complex ip, scale;
    int i /*, j */;
 
    if ( M==ZMNULL || hh==ZVNULL )
        error(E_NULL,"zhhtrrows");
    if ( M->n != hh->dim )
        error(E_RANGE,"zhhtrrows");
    if ( i0 < 0 || i0 > M->m || j0 < 0 || j0 > M->n )
        error(E_BOUNDS,"zhhtrrows");
 
    if ( beta == 0.0 )  return (M);
 
    /* for each row ... */
    for ( i = i0; i < M->m; i++ )
    {   /* compute inner product */
        ip = __zip__(&(M->me[i][j0]),&(hh->ve[j0]),
                 (int)(M->n-j0),Z_NOCONJ);
        /**************************************************
        ip = 0.0;
        for ( j = j0; j < M->n; j++ )
            ip += M->me[i][j]*hh->ve[j];
        **************************************************/
        scale.re = -beta*ip.re;
        scale.im = -beta*ip.im;
        /* if ( scale == 0.0 ) */
        if ( is_zero(scale) )
            continue;
 
        /* do operation */
        __zmltadd__(&(M->me[i][j0]),&(hh->ve[j0]),scale,
                (int)(M->n-j0),Z_CONJ);
        /**************************************************
        for ( j = j0; j < M->n; j++ )
            M->me[i][j] -= scale*hh->ve[j];
        **************************************************/
    }
 
    return (M);
}
 
 
/* zhhtrcols -- transform a matrix by a Householder vector by columns
    starting at row i0 from column j0 -- in-situ */
ZMAT    *zhhtrcols(M,i0,j0,hh,beta)
ZMAT    *M;
int i0, j0;
ZVEC    *hh;
double  beta;
{
    /* Real ip, scale; */
    complex scale;
    int i /*, k */;
    static  ZVEC    *w = ZVNULL;
 
    if ( M==ZMNULL || hh==ZVNULL )
        error(E_NULL,"zhhtrcols");
    if ( M->m != hh->dim )
        error(E_SIZES,"zhhtrcols");
    if ( i0 < 0 || i0 > M->m || j0 < 0 || j0 > M->n )
        error(E_BOUNDS,"zhhtrcols");
 
    if ( beta == 0.0 )  return (M);
 
    w = zv_resize(w,M->n);
    MEM_STAT_REG(w,TYPE_ZVEC);
    zv_zero(w);
 
    for ( i = i0; i < M->m; i++ )
        /* if ( hh->ve[i] != 0.0 ) */
        if ( ! is_zero(hh->ve[i]) )
        __zmltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
                (int)(M->n-j0),Z_CONJ);
    for ( i = i0; i < M->m; i++ )
        /* if ( hh->ve[i] != 0.0 ) */
        if ( ! is_zero(hh->ve[i]) )
        {
        scale.re = -beta*hh->ve[i].re;
        scale.im = -beta*hh->ve[i].im;
        __zmltadd__(&(M->me[i][j0]),&(w->ve[j0]),scale,
                (int)(M->n-j0),Z_CONJ);
        }
    return (M);
}