ltfat_lapack.c

Matlab时频分析工具箱,希望能对大家有所帮助啊

#include "config.h"
#include "cblas.h"

/* Use the LAPACK library supplied with ATLAS */
#ifdef HAVE_ATLASLAPACK
#include "clapack.h"
#endif /* end of HAVE_ATLASLAPACK */

/* Call Fortran LAPACK */
#ifdef HAVE_LAPACK
#include "f77-fcn.h"


#ifdef __cplusplus
extern "C" {
#endif

  F77_RET_T
  F77_FUNC (zposv, ZPOSV) (F77_CONST_CHAR_ARG_DECL uplo,
			   const int *n, const int *nrhs, 
			   double *a, const int *lda,
			   double *b, const int *ldb, 
			   int *info
			   F77_CHAR_ARG_LEN_DECL);


  F77_RET_T
  F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL jobu,
                             F77_CONST_CHAR_ARG_DECL jobvt,
                             const int *M, const int *N,
			     double* A, const int *lda,
			     double *S, double* U, const int *ldu,
                             double *VT, const int *ldvt, double *work,
			     const int *lwork,
                             double *rwork, int *info
                             F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);
#ifdef __cplusplus
}
#endif

#endif /* end of HAVE_LAPACK */

/* ----- Compute Cholesky factorization  ------------
 *
 * For simplification, the interface assumes that we
 * are using column-major format and storing the upper
 * triangle
 */
int ltfat_zposv(const int N, const int NRHS,
		const ltfat_complex *A, const int lda,
		ltfat_complex *B, const int ldb)
#ifdef HAVE_CBLASLAPACK

{
  return clapack_zposv(CblasColMajor, CblasUpper, N, NRHS, (double*)A, lda,
		       (double*)B, ldb);
}
#endif
#ifdef HAVE_LAPACK
{
  int info;
  char u;

  u = 'U';

  F77_FUNC (zposv, ZPOSV) (F77_CONST_CHAR_ARG2 (&u, 1),
			   &N, &NRHS, (double*)A, &lda,
			   (double*)B, &ldb,
			   &info
			   F77_CHAR_ARG_LEN (1)
			   );

  return info;
}
#endif

/* ----- Compute SVD factorization  ------------ */
int ltfat_zgesvd(const int M, const int N,
		 ltfat_complex *A, const int lda,
		 double *S, ltfat_complex *U, const int ldu,
		 ltfat_complex *VT, const int ldvt)
#ifdef HAVE_LAPACK
{

  int lwork, info, minMN, maxMN;
  double workquery[2]; 
  double *rwork, *work;

  char jobu;
  char jobvt;

  /* Set constants to declare thin SVD */
  jobu = 'S';
  jobvt = 'S';  

  minMN = M < N ? M : N;
  maxMN = M > N ? M : N;

  /* Allocate workspace */
  rwork = (double*)ltfat_malloc(5*maxMN*sizeof(double));

  /* Ask ZGESVD what the dimension of WORK should be. */
  lwork = -1;
  F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG2 (&jobu, 1),
                             F77_CONST_CHAR_ARG2 (&jobvt, 1),
                             &M, &N, (double*)A, &lda, S, (double*)U,
			     &ldu, (double*)VT,
                             &ldvt, (double*)(&workquery), &lwork,
                             rwork, &info
                             F77_CHAR_ARG_LEN (1)
                             F77_CHAR_ARG_LEN (1));

  /* Get the result from real part of work */
  lwork = (int)(workquery[0]);

  /* Allocate more workspace */
  work = (double*)ltfat_malloc(lwork*sizeof(ltfat_complex));
  
  /* Call the function for real this time */

  F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG2 (&jobu, 1),
                             F77_CONST_CHAR_ARG2 (&jobvt, 1),
                             &M, &N, (double*)A, &lda, S,
			     (double*)U, &ldu, (double*)VT,
                             &ldvt, work, &lwork,
                             rwork, &info
                             F77_CHAR_ARG_LEN (1)
                             F77_CHAR_ARG_LEN (1));


  /* Free workspace */
  ltfat_free(rwork);
  ltfat_free(work);

  return info;

}
#endif