hholder.c File Reference

Implementation and use of Householder transform. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "tpcextensions.h"
#include "tpclinopt.h"

Go to the source code of this file.

Functions
double	householder_transform (double *v, int N)
int	householder_hm (double tau, double *vector, double **matrix, int rowNr, int columnNr)
int	householder_hv (double tau, int size, double *v, double *w)
double	householder_norm (double *v, int size)

Detailed Description

Implementation and use of Householder transform.

Author

Kaisa Liukko, Vesa Oikonen

These routines are based on the code provided in the GSL library (https://sources.redhat.com/gsl/).

Todo

Use const variables when possible.

Definition in file hholder.c.

Function Documentation

householder_hm()

int householder_hm	(	double	tau,
		double *	vector,
		double **	matrix,
		int	rowNr,
		int	columnNr )

Applies a householder transformation defined by vector "vector" and scalar tau to the left-hand side of the matrix. (I - tau vector vector^T)*matrix The result of the transform is stored in matrix.

Returns

Returns 0 if ok.

Parameters

tau	Coefficient defining householder transform.
vector	Vector defining householder transform (of size rowNr).
matrix	the matrix that is to be transformed.
rowNr	Nr of rows in matrix.
columnNr	Nr of columns in matrix.

Definition at line 87 of file hholder.c.

  {
  int i, j;
  double wj;
 
  if(tau==0.0) return(0); // success
  if(rowNr<1 || columnNr<1) return(1);
  for(j=0; j<columnNr; j++) {
    /* Compute wj = vk Akj */
    wj=matrix[0][j];
    for(i=1; i<rowNr; i++)  /* note, computed for v(0) = 1 above */
      wj += vector[i]*matrix[i][j];
    /* Aij = Aij - tau vi wj */
    /* i = 0 */
    matrix[0][j]-=tau*wj;
    /* i = 1 .. M-1 */
    for(i=1; i<rowNr; i++) matrix[i][j]-=tau*vector[i]*wj;
  }
  return(0);
}

Referenced by qr_decomp().

householder_hv()

int householder_hv	(	double	tau,
		int	size,
		double *	v,
		double *	w )

Applies a householder transformation defined by vector v and coefficient tau to vector w w = (I - tau v v^T) w.

Returns

Returns 0 if ok.

Parameters

tau	Coefficient defining householder transform.
size	Size of vectors v and w.
v	Vector v.
w	Vector w.

Definition at line 126 of file hholder.c.

  {
  int i;
  double d;
 
  if(tau==0) return(0); // success
  if(size<1) return(1);
  /* d = v'w */
  d=w[0]; for(i=1; i<size; i++) d+=v[i]*w[i];
  /* w = w - tau (v) (v'w) */
  w[0]-=tau*d;
  for(i=1; i<size; i++) w[i]-=tau*v[i]*d;
  return(0);
}

householder_norm()

double householder_norm	(	double *	v,
		int	size )

Calculates the euclidean norm of vector v[].

Returns

Returns the euclidean norm of a vector.

Parameters

v	Vector v.
size	Size of vector v[].

Definition at line 155 of file hholder.c.

  {
  double ss=0.0;
  int i;
 
  for(i=0; i<size; i++) ss+=v[i]*v[i];
  return sqrt(ss);
}

householder_transform()

double householder_transform	(	double *	v,
		int	N )

This function prepares a Householder transformation P = I - tau h h^T which can be used to zero all the elements of the input vector except the first one that will get value beta. On output the elements 1 - size-1 of the vector h are stored in locations vector[1] - vector[size-1] of the input vector and value of beta is stored in location vector[0].

Returns

The scalar tau is returned.

Parameters

v	The N-vector to be transformed.
N	size of the vector.

Definition at line 33 of file hholder.c.

  {
  double vnorm, alpha, beta, tau;
  int n;
 
  if(N<1) return 0.0; // tau = 0
  /* Euclidean norm of the vector starting from the second value */
  vnorm=0.0; for(n=1; n<N; n++) vnorm+=v[n]*v[n];
  vnorm=sqrt(vnorm); if(isnan(vnorm) || vnorm==0.0) return 0.0; // tau = 0
 
  /* Computing the coefficient tau */
  alpha=v[0];
  beta= - (alpha >= 0.0 ? +1.0 : -1.0) * hypot(alpha, vnorm);
  tau=(beta-alpha)/beta ;
 
  /* Scale the Householder vector so that the first element will be 1.
   * (Scaling is also affecting the coefficient tau).
   * Without scaling, the first element would have value (alpha - beta). */
#if(0) // Kaisa's version
  {
    double os=1.0/(alpha-beta);
    for(n=1; n<N; n++) v[n]*=os;
    v[0]=beta;
  }
#else // Vesa's version
  {
    double s=alpha-beta;
    if(fabs(s)>DBL_MIN) {
      v[0]=beta;
      for(n=1; n<N; n++) v[n]*=(1.0/s);
    } else {
      v[0]=beta;
      for(n=1; n<N; n++) v[n]*=(doubleMachEps()/s);
      for(n=1; n<N; n++) v[n]*=(1.0/doubleMachEps());
    }
  }
#endif
 
  return tau;
}

Referenced by qr_decomp().