home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: 22 gnu
/
22-gnu.zip
/
mesch12a.zip
/
zhessen.c
< prev
next >
Wrap
C/C++ Source or Header
|
1994-01-13
|
4KB
|
153 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.
**
***************************************************************************/
/*
File containing routines for determining Hessenberg
factorisations.
Complex version
*/
static char rcsid[] = "$Id: zhessen.c,v 1.1 1994/01/13 04:27:00 des Exp $";
#include <stdio.h>
#include "zmatrix.h"
#include "zmatrix2.h"
/* zHfactor -- compute Hessenberg factorisation in compact form.
-- factorisation performed in situ
-- for details of the compact form see zQRfactor.c and zmatrix2.doc */
ZMAT *zHfactor(A, diag)
ZMAT *A;
ZVEC *diag;
{
static ZVEC *tmp1 = ZVNULL;
Real beta;
int k, limit;
if ( ! A || ! diag )
error(E_NULL,"zHfactor");
if ( diag->dim < A->m - 1 )
error(E_SIZES,"zHfactor");
if ( A->m != A->n )
error(E_SQUARE,"zHfactor");
limit = A->m - 1;
tmp1 = zv_resize(tmp1,A->m);
MEM_STAT_REG(tmp1,TYPE_ZVEC);
for ( k = 0; k < limit; k++ )
{
zget_col(A,k,tmp1);
zhhvec(tmp1,k+1,&beta,tmp1,&A->me[k+1][k]);
diag->ve[k] = tmp1->ve[k+1];
/* printf("zHfactor: k = %d, beta = %g, tmp1 =\n",k,beta);
zv_output(tmp1); */
zhhtrcols(A,k+1,k+1,tmp1,beta);
zhhtrrows(A,0 ,k+1,tmp1,beta);
/* printf("# at stage k = %d, A =\n",k); zm_output(A); */
}
return (A);
}
/* zHQunpack -- unpack the compact representation of H and Q of a
Hessenberg factorisation
-- if either H or Q is NULL, then it is not unpacked
-- it can be in situ with HQ == H
-- returns HQ
*/
ZMAT *zHQunpack(HQ,diag,Q,H)
ZMAT *HQ, *Q, *H;
ZVEC *diag;
{
int i, j, limit;
Real beta, r_ii, tmp_val;
static ZVEC *tmp1 = ZVNULL, *tmp2 = ZVNULL;
if ( HQ==ZMNULL || diag==ZVNULL )
error(E_NULL,"zHQunpack");
if ( HQ == Q || H == Q )
error(E_INSITU,"zHQunpack");
limit = HQ->m - 1;
if ( diag->dim < limit )
error(E_SIZES,"zHQunpack");
if ( HQ->m != HQ->n )
error(E_SQUARE,"zHQunpack");
if ( Q != ZMNULL )
{
Q = zm_resize(Q,HQ->m,HQ->m);
tmp1 = zv_resize(tmp1,H->m);
tmp2 = zv_resize(tmp2,H->m);
MEM_STAT_REG(tmp1,TYPE_ZVEC);
MEM_STAT_REG(tmp2,TYPE_ZVEC);
for ( i = 0; i < H->m; i++ )
{
/* tmp1 = i'th basis vector */
for ( j = 0; j < H->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(HQ,j,tmp2);
r_ii = zabs(tmp2->ve[j+1]);
tmp2->ve[j+1] = diag->ve[j];
tmp_val = (r_ii*zabs(diag->ve[j]));
beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
/* printf("zHQunpack: j = %d, beta = %g, tmp2 =\n",
j,beta);
zv_output(tmp2); */
zhhtrvec(tmp2,beta,j+1,tmp1,tmp1);
}
/* insert into Q */
zset_col(Q,i,tmp1);
}
}
if ( H != ZMNULL )
{
H = zm_copy(HQ,H);
limit = H->m;
for ( i = 1; i < limit; i++ )
for ( j = 0; j < i-1; j++ )
H->me[i][j].re = H->me[i][j].im = 0.0;
}
return HQ;
}