hholder.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
/*****************************************************************************/
#include "tpclinopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
double householder_transform(
  double *v,
  int N
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/
int householder_hm (
  double tau,
  double *vector,
  double **matrix,
  int rowNr,
  int columnNr
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int householder_hv(
  double tau,
  int size,
  double *v,
  double *w
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
double householder_norm(
  double *v,
  int size
) {
  double ss=0.0;
  int i;
 
  for(i=0; i<size; i++) ss+=v[i]*v[i];
  return sqrt(ss);
}
/*****************************************************************************/
 
/*****************************************************************************/