appl.c

计算稀疏正交因子的matlab算法库

/* Apply Householder vectors obtained from Householder QR.
   Mikael Adlers, 
   Linkoping University, Sweden
   [e-mail]
   990118

   This routine should not be used seperately. It is part of the sqr2 package.
   If you encounter any problems when compiling this file, please contact the
   author.

   Usage:
   B=appL(b,H,tau,mode);
   
   b : the input vector to be transformed
   H : lower trapezoidal part contains Householder transformation.
       The transformation is scaled so first element is equal to one
   tau : scale factor 
   mode : 1 = notranspose (H*b). 2 = transpose (H'*b)

   B=(I-tau_i v_i v_i^T)*b  for all i= 1, .., n
   v_i is the columns in H + diagonal equal to 1.

*/
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "matrix.h"
#include "mex.h"

#define A_in       prhs[0]
#define H_in       prhs[1]
#define tau_in     prhs[2]
#define mode_in    prhs[3]

#define B_out      plhs[0]

/* ------------------- Entry point ------------------------ */

void mexFunction
(
/* === Parameters ======================================================= */

    int nlhs,			/* number of left-hand sides */
    mxArray *plhs [],		/* left-hand side matrices */
    int nrhs,			/* number of right--hand sides */
    const mxArray *prhs []	/* right-hand side matrices */
)
{
  unsigned int m1,n1,nH,mH;
  mxArray *y;
  double *tau;
  double *work,*H,*B,*A,*col,w;
  int *H_Ir,*H_Jc;
  int lwork;
  int i,j,k,info,lda,mode;
  char s1,s2;
  
  /* --- Check that the input paramenters are correct --- */      
  
  if ((nrhs > 5) || (nrhs<4)) 
    {mexErrMsgTxt("Four input arguments required.");}
  if (nlhs != 1) 
    {mexErrMsgTxt("One output argument required");}
  
  m1=mxGetM(A_in);
  n1=mxGetN(A_in);
  mH=mxGetM(H_in);
  nH=mxGetN(H_in);
  
  if(!mxIsDouble(A_in) || mxIsComplex(A_in))
    {mexErrMsgTxt("First argument must be must be a real numeric matrix");}
    
  if ((m1>0)&&(n1>0)){
      if (!mxIsSparse(H_in))
      	{
	 mexErrMsgTxt("Only sparse Householdervectors allowed.");
      	}
      else
	{
	  B_out   = mxCreateDoubleMatrix(m1,n1,mxREAL);

	  H   = (double *)mxGetPr(H_in);
	  H_Ir= (int *)   mxGetIr(H_in);
	  H_Jc= (int *)   mxGetJc(H_in);
	  A   = (double *)mxGetPr(A_in);
	  tau = (double *)mxGetPr(tau_in);
	  B   = (double *)mxGetPr(B_out);
	  mode= (int) mxGetScalar(mode_in);

	  for (i=0;i<m1*n1;i++)
	      B[i]=A[i];

	  /*  dgeqrf_(&m1,&n1,R,&lda,tau,work,&lwork,&info);*/
	  if (mode==2){ /* Transpose mode H'*b */
	    /*dormqr_("L","N",&m1,&n1,&nH,H,&mH,tau,B,&m1,work,&lwork,&info);*/
	    /* nH is the number of transformations */
	    for (i=0;i<nH;i++){
	      
	      /* Transformation i affects only b(i:m1,:) */
	      for (j=0;j<n1;j++){
		col=B+(j*m1); /* pointer to b(:,j) */

		
		w=col[i]; 
		for (k=H_Jc[i];k<H_Jc[i+1];k++) w+=H[k]*col[H_Ir[k]];

		/*for (k=i+1;k<m;k++) w+=v[k]*col[k];*/
     		if (w!=0) {
		  
		  w*=tau[i];
		  col[i]-=w;
		  for(k=H_Jc[i];k<H_Jc[i+1];k++) col[H_Ir[k]]-=H[k]*w;
		}
		
	      }
	    }
	  }
	  else { /* No transpose mode H*b */
	  /*dormqr_("L","T",&m1,&n1,&nH,H,&mH,tau,B,&m1,work,&lwork,&info);*/
	    for (i=nH-1;i>=0;i--){
	      
	      /* Transformation i affects only b(i:m1,:) */
	      for (j=0;j<n1;j++){
		col=B+(j*m1); /* pointer to b(:,j) */
		
		
		w=col[i]; 
		for (k=H_Jc[i];k<H_Jc[i+1];k++) w+=H[k]*col[H_Ir[k]];
		
		/*for (k=i+1;k<m;k++) w+=v[k]*col[k];*/
     		if (w!=0) {
		  
		  w*=tau[i];
		  col[i]-=w;
		  for(k=H_Jc[i];k<H_Jc[i+1];k++) col[H_Ir[k]]-=H[k]*w;
		}
		
	      }
	      
	    }
	    
	  }
	}
  }
  else
    {
      B_out   = mxCreateDoubleMatrix(m1,n1,mxREAL);
    }
}