#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "liblapackwrap/lapackwrap_real_double.h"
#include "liblapackwrap/lapackwrap_functions_real_double.c"

Include dependency graph for lapackwrap_real_double.c:

Defines
#define	ACML 0
#define	CELL 1
#define	FCALL
Functions
int	llw_real_double_stegr_a (char jobz, char range, int N, LLW_CONST LLW_real_double_REAL D, LLW_CONST LLW_real_double_REAL E, LLW_real_double_REAL VL, LLW_real_double_REAL VU, int IL, int IU, int M, LLW_real_double_REAL W, LLW_real_double_COMPLEX Z, int LDZ, int **ISUPPZ)
int	llw_real_double_stegr (char jobz, char range, int N, LLW_CONST LLW_real_double_REAL D, LLW_CONST LLW_real_double_REAL E, LLW_real_double_REAL VL, LLW_real_double_REAL VU, int IL, int IU, int M, LLW_real_double_REAL W, LLW_real_double_COMPLEX Z, int LDZ, int ISUPPZ)
int	llw_real_double_sbtrd_a (char vect, char uplo, int N, int KD, LLW_CONST LLW_real_double_COMPLEX A, int LDAB, LLW_real_double_REAL D, LLW_real_double_REAL E, LLW_real_double_COMPLEX *Q)
int	llw_real_double_sbtrd (char vect, char uplo, int N, int KD, LLW_CONST LLW_real_double_COMPLEX AB, int LDAB, LLW_real_double_REAL D, LLW_real_double_REAL E, LLW_real_double_COMPLEX Q, int LDQ)
int	llw_real_double_sytrd (char uplo, int N, LLW_CONST LLW_real_double_COMPLEX A, LLW_real_double_REAL D, LLW_real_double_REAL E, LLW_real_double_COMPLEX TAU, LLW_real_double_COMPLEX **a_out)
int	llw_real_double_orgtr (char uplo, int N, LLW_real_double_COMPLEX A, LLW_real_double_COMPLEX TAU)

Define Documentation

#define ACML 0

Definition at line 27 of file lapackwrap_real_double.c.

#define CELL 1

Definition at line 28 of file lapackwrap_real_double.c.

#define FCALL

Definition at line 32 of file lapackwrap_real_double.c.

Function Documentation

int llw_real_double_orgtr	(	char	uplo,
		int	N,
		LLW_real_double_COMPLEX *	A,
		LLW_real_double_COMPLEX *	TAU
	)

Turns the output of a tridiagonalization into a matrix which will help convert the eigvals to the old basis

Returns:: INFO, 100*the following error values: 1: Failed to allocate space.

Note:: A is NOT LLW_CONST. orgtr will overwrite it. This is probably what you want, so it's preserved.

Definition at line 270 of file lapackwrap_real_double.c.

References LAPACK_COMPLEX_TYPENAME, LLW_real_double_COMPLEX, and ORGTR.

                                                                                                      {
  int info;
  #if LAPACK == ACML
  ORGTR(uplo, N, (LAPACK_COMPLEX_TYPENAME*)A, N, (LAPACK_COMPLEX_TYPENAME*)TAU, &info);
  #else
  LLW_real_double_COMPLEX opt_work;
  #ifdef FCALL
  int MONE = -1;
  ORGTR(&uplo, &N,  (LAPACK_COMPLEX_TYPENAME*)A, &N,  (LAPACK_COMPLEX_TYPENAME*)TAU,  (LAPACK_COMPLEX_TYPENAME*)&opt_work, &MONE, &info);
  #else
  ORGTR(uplo, N,  (LAPACK_COMPLEX_TYPENAME*)A, N,  (LAPACK_COMPLEX_TYPENAME*)TAU,  (LAPACK_COMPLEX_TYPENAME*)&opt_work, -1, &info);
  #endif
  if(info != 0) {
    return info;
  }
  int OWORK = (int)
      #ifdef LLW_real_double_IS_COMPLEX
              opt_work.re
      #else
              opt_work
            #endif
    ;
  LLW_real_double_COMPLEX* work = (LLW_real_double_COMPLEX*)malloc(sizeof(LLW_real_double_COMPLEX)*OWORK);
  #ifdef FCALL
  ORGTR(&uplo, &N,  (LAPACK_COMPLEX_TYPENAME*)A, &N,  (LAPACK_COMPLEX_TYPENAME*)TAU,  (LAPACK_COMPLEX_TYPENAME*)work, &OWORK, &info);
  #else
  ORGTR(uplo, N,  (LAPACK_COMPLEX_TYPENAME*)A, N,  (LAPACK_COMPLEX_TYPENAME*)TAU,  (LAPACK_COMPLEX_TYPENAME*)work, OWORK, &info);
  #endif
  #endif
  return info;
}

int llw_real_double_sbtrd	(	char	vect,
		char	uplo,
		int	N,
		int	KD,
		LLW_CONST LLW_real_double_COMPLEX *	AB,
		int	LDAB,
		LLW_real_double_REAL *	D,
		LLW_real_double_REAL *	E,
		LLW_real_double_COMPLEX *	Q,
		int	LDQ
	)

Gets the tri-diagonal matrix for a banded matrix.

Definition at line 188 of file lapackwrap_real_double.c.

References LAPACK_COMPLEX_TYPENAME, LLW_real_double_COMPLEX, and SBTRD.

Referenced by llw_real_double_sbtrd_a().

                                                                                                                                                                                                       {
  LLW_real_double_COMPLEX *ab;
  #ifdef LLW_CONST
  ab = (LLW_real_double_COMPLEX*)malloc(N*(KD+1)*sizeof(LLW_real_double_COMPLEX));
  if(ab == NULL) return 100;
  memcpy(ab, AB, N*(KD+1)*sizeof(LLW_real_double_COMPLEX));
  #else
  ab = AB;
  #endif
  int info;
  #if LAPACK == ACML
  SBTRD(vect, uplo, N, KD, (LAPACK_COMPLEX_TYPENAME*)ab, LDAB, D, E,  (LAPACK_COMPLEX_TYPENAME*)Q, LDQ, &info);
  #else
  LLW_real_double_COMPLEX *work = (LLW_real_double_COMPLEX*)malloc(N*sizeof(LLW_real_double_COMPLEX));
  if(work == NULL) {
    #ifdef LLW_CONST
    free(ab);
    #endif
    return 100;
  }
  #ifdef FCALL
  SBTRD(&vect, &uplo, &N, &KD,  (LAPACK_COMPLEX_TYPENAME*)ab, &LDAB, D, E,  (LAPACK_COMPLEX_TYPENAME*)Q, &LDQ,  (LAPACK_COMPLEX_TYPENAME*)work, &info);
  #else
  SBTRD(vect, uplo, N, KD,  (LAPACK_COMPLEX_TYPENAME*)ab, LDAB, D, E,  (LAPACK_COMPLEX_TYPENAME*)Q, LDQ,  (LAPACK_COMPLEX_TYPENAME*)work, &info);
  #endif
  free(work);
  #endif
  return info;
}

Here is the caller graph for this function:

int llw_real_double_sbtrd_a	(	char	vect,
		char	uplo,
		int	N,
		int	KD,
		LLW_CONST LLW_real_double_COMPLEX *	A,
		int	LDAB,
		LLW_real_double_REAL **	D,
		LLW_real_double_REAL **	E,
		LLW_real_double_COMPLEX **	Q
	)

This is the same as llw_real_double_sbtrd, but with dynamic allocation of D, E, and Q.

Definition at line 171 of file lapackwrap_real_double.c.

References LLW_real_double_COMPLEX, LLW_real_double_REAL, and llw_real_double_sbtrd().

                                                                                                                                                                                                  {
  *D = (LLW_real_double_REAL*)malloc(N*sizeof(LLW_real_double_REAL));
  if(*D == NULL) return 100;
  *E = (LLW_real_double_REAL*)malloc((N-1)*sizeof(LLW_real_double_REAL));
  if(*E == NULL) {
    free(D);
    return 100;
  }
  *Q = (LLW_real_double_COMPLEX*)malloc(N*N*sizeof(LLW_real_double_COMPLEX));
  if(*Q == NULL) {
    free(D);
    free(E);
    return 100;
  }
  return llw_real_double_sbtrd(vect, uplo, N, KD, A, LDAB, *D, *E, *Q, N);
}

Here is the call graph for this function:

int llw_real_double_stegr	(	char	jobz,
		char	range,
		int	N,
		LLW_CONST LLW_real_double_REAL *	D,
		LLW_CONST LLW_real_double_REAL *	E,
		LLW_real_double_REAL	VL,
		LLW_real_double_REAL	VU,
		int	IL,
		int	IU,
		int *	M,
		LLW_real_double_REAL *	W,
		LLW_real_double_COMPLEX *	Z,
		int	LDZ,
		int *	ISUPPZ
	)

wrapper around dstegr/dstegr_

jobz passed straight into LAPACK
range passed straight into LAPACK
N passed straight into LAPACK
D passed straight into LAPACK
E passed straight into LAPACK
VL passed straight into LAPACK
VU passed straight into LAPACK
IL passed straight into LAPACK
IU passed straight into LAPACK
M passed straight into LAPACK
W passed straight into LAPACK
Z passed straight into LAPACK
LDZ passed straight into LAPACK
ISUPPZ passed straight into LAPACK
Returns:
INFO, 100*the following error values: 1: Failed to allocate space.

Definition at line 68 of file lapackwrap_real_double.c.

References LAPACK_COMPLEX_TYPENAME, LLW_real_double_REAL, and STEGR.

Referenced by llw_real_double_stegr_a().

                                                                                                                                                                                                                                                                                   {
  /*
  fprintf(stderr, "llw_real_double_stegr: %c %c %d %p %p %g %g %d %d %p %p %p %d %p\n", jobz, range, N, (void*)D, (void*)E, VL, VU, IL, IU, (void*)M, (void*)W, (void*)Z, LDZ, (void*)ISUPPZ);
  */
  /*Now do our STEGR call*/
  /*NOTE e_in is of length N. We copy only N-1 (since E is N-1 values; it's the first subdiagonal. HOWEVER, e_in[N-1] (i.e., the Nth element) is used internally by LAPACK as a workspace, so we need it in here.*/
  LLW_real_double_REAL *e_in;
  LLW_real_double_REAL *d_in;
  #ifdef LLW_CONST
  e_in = malloc(N*sizeof(LLW_real_double_REAL));
  d_in = malloc(N*sizeof(LLW_real_double_REAL));
  if((e_in == NULL) || (d_in == NULL)) {
    if(e_in != NULL) free(e_in);
    if(d_in != NULL) free(d_in);
    return 100;
  }
  memcpy(e_in, E, (N-1)*sizeof(LLW_real_double_REAL));
  memcpy(d_in, D, N*sizeof(LLW_real_double_REAL));
  #else
  e_in = E;
  d_in = D;
  #endif
  /*Initialize our output storage*/
  /*
  *M = 0;
  for(int i=0; i<N; i++) {
    W[i] = 0;
    for(int j=0; j<LDZ; j++) {
      Z[j + i*LDZ] = 0;
    }
    ISUPPZ[2*i] = ISUPPZ[2*i + 1] = 0;
  }
  for(int i=0; i<N-1; i++) {
    fprintf(stderr, "d[%d]=", i);
    fprintf(stderr, "%g ", d_in[i]);
    fprintf(stderr, "e[%d]=", i);
    fprintf(stderr, "=%g\n", e_in[i]);
    fprintf(stderr, "D[%d]=", i);
    fprintf(stderr, "%g ", D[i]);
    fprintf(stderr, "E[%d]=", i);
    fprintf(stderr, "=%g\n", E[i]);
  }
  fprintf(stderr, "D[%d]=%g\n", N, D[N-1]);
  fprintf(stderr, "d[%d]=%g\n", N, d_in[N-1]);
  */
  /*First, query optimal work size*/
  /*
  for(int i=0; i<N; i++) {
    fprintf(stderr, "D[%d]=%g d_in=%g E=%g e_in=%g\n", 
      i, 
      D[i], 
      d_in[i], 
      E[i], 
      (i == (N-1))?0:e_in[1]);
  }
  */
  int info;
#if LAPACK == ACML
  /*ACML's C interface doesn't require the work and iwork arrays.
   *So with ACML, we just dive right in.
   */
  STEGR(jobz, range, N, d_in, e_in, VL, VU, IL, IU, 1e-12, M, W, (LAPACK_COMPLEX_TYPENAME*)Z, LDZ, ISUPPZ, &info);
#else
  /*The rest of the LAPACK implementations require work and iwork,
   *  so we'll query the optimal values before making the real call.
   */
  LLW_real_double_REAL opt_lwork;
  int opt_liwork;
#ifdef FCALL
  /*Look at all the fun in a FORTRAN call. Whee!*/
  LLW_real_double_REAL EPSILON=1e-12;
  int MONE = -1;
  STEGR(&jobz, &range, &N, d_in, e_in, &VL, &VU, &IL, &IU, &EPSILON, M, W, (LAPACK_COMPLEX_TYPENAME*)Z, &LDZ, ISUPPZ, &opt_lwork, &MONE, &opt_liwork, &MONE, &info);
#else
  STEGR(jobz, range, N, d_in, e_in, VL, VU, IL, IU, 1e-12, M, W, (LAPACK_COMPLEX_TYPENAME*)Z, LDZ, ISUPPZ, &opt_lwork, -1, &opt_liwork, -1, &info);
#endif
  if(info != 0) {
    #ifdef LLW_CONST
    free(e_in);
    free(d_in);
    #endif
    return info;
  }
  LLW_real_double_REAL *work = malloc(sizeof(LLW_real_double_REAL)*((int)opt_lwork));
  int *iwork = malloc(sizeof(int)*opt_liwork);
  #ifdef FCALL
  int LWORK = opt_lwork;
  STEGR(&jobz, &range, &N, d_in, e_in, &VL, &VU, &IL, &IU, &EPSILON, M, W, (LAPACK_COMPLEX_TYPENAME*)Z, &LDZ, ISUPPZ, work, &LWORK, iwork, &opt_liwork, &info);
  #else
  STEGR(jobz, range, N, d_in, e_in, VL, VU, IL, IU, 1e-12, M, W, (LAPACK_COMPLEX_TYPENAME*)Z, LDZ, ISUPPZ, work, opt_lwork, iwork, opt_liwork, &info);
  #endif
  free(work);
  free(iwork);
  #endif
  #ifdef LLW_CONST
  free(e_in);
  free(d_in);
  #endif
  return info;
}

Here is the caller graph for this function:

int llw_real_double_stegr_a	(	char	jobz,
		char	range,
		int	N,
		LLW_CONST LLW_real_double_REAL *	D,
		LLW_CONST LLW_real_double_REAL *	E,
		LLW_real_double_REAL	VL,
		LLW_real_double_REAL	VU,
		int	IL,
		int	IU,
		int *	M,
		LLW_real_double_REAL **	W,
		LLW_real_double_COMPLEX **	Z,
		int *	LDZ,
		int **	ISUPPZ
	)

same as dstegr, but allocation of W, Z, and SUPPZ is done internally

jobz passed straight into LAPACK
range passed straight into LAPACK
N passed straight into LAPACK
D passed straight into LAPACK
E passed straight into LAPACK
VL passed straight into LAPACK
VU passed straight into LAPACK
IL passed straight into LAPACK
IU passed straight into LAPACK
M passed straight into LAPACK
W pointer to a pointer to a LLW_real_double_REAL; will be used to store eigvals
Z pointer to a pointer to a LLW_real_double_REAL; will be used to store eigvecs (column-major)
LDZ pointer to an int; will be used to store number of eigvecs in Z
ISUPPZ pointer to pointer to a LLW_real_double_REAL; will be used to store support of matrix.
Returns:
INFO, 100*the following error values: 1: Failed to allocate space.

Definition at line 35 of file lapackwrap_real_double.c.

References LLW_real_double_COMPLEX, LLW_real_double_REAL, and llw_real_double_stegr().

                                                                                                                                                                                                                                                                                         {
  *W = malloc(sizeof(LLW_real_double_REAL)*N);
  if(*W == NULL) {
    return 100;
  }
  /*Right now, we're just doing N.*/
  *LDZ = N;
  *Z = malloc(sizeof(LLW_real_double_COMPLEX)*(*LDZ)*N);
  if(*Z == NULL) {
    free(*W);
    *W = NULL;
    return 100;
  }
  *ISUPPZ = malloc(sizeof(LLW_real_double_REAL)*2*N);
  if(*ISUPPZ == NULL) {
    free(*W);
    *W = NULL;
    free(*Z);
    *Z = NULL;
    return 100;
  }
  int info = llw_real_double_stegr(jobz, range, N, D, E, VL, VU, IL, IU, M, *W, *Z, *LDZ, *ISUPPZ);
  if(info != 0) {
    free(*W);
    *W = NULL;
    free(*Z);
    *Z = NULL;
    free(*ISUPPZ);
    *ISUPPZ = NULL;
  }
  return info;
}

Here is the call graph for this function:

int llw_real_double_sytrd	(	char	uplo,
		int	N,
		LLW_CONST LLW_real_double_COMPLEX *	A,
		LLW_real_double_REAL *	D,
		LLW_real_double_REAL *	E,
		LLW_real_double_COMPLEX *	TAU,
		LLW_real_double_COMPLEX **	a_out
	)

LAPACK wrapper for DSYTRD.

Returns:: the resulting INFO

Definition at line 223 of file lapackwrap_real_double.c.

References LAPACK_COMPLEX_TYPENAME, LLW_real_double_COMPLEX, and SYTRD.

                                                                                                                                                                                                   {
  LLW_real_double_COMPLEX **aa;
  LLW_real_double_COMPLEX *a;
  if(a_out != NULL) {
    a = *a_out;
  }else{
    /*Point aa to a LLW_real_double_COMPLEX pointer on the stack if the user doesn't want it*/
    aa = &a;
  }
  *aa = malloc(sizeof(LLW_real_double_COMPLEX)*N*N);
  if(a == NULL) return 100;
  memcpy(*aa, A, N*N*sizeof(LLW_real_double_COMPLEX));
  int info;
  #if LAPACK == ACML
  SYTRD(uplo, N,  (LAPACK_COMPLEX_TYPENAME*)*aa, N, D, E,  (LAPACK_COMPLEX_TYPENAME*)TAU, &info);
  #else
  LLW_real_double_COMPLEX opt_work;
  #ifdef FCALL
  int MONE = -1;
  SYTRD(&uplo, &N,  (LAPACK_COMPLEX_TYPENAME*)*aa, &N, D, E,  (LAPACK_COMPLEX_TYPENAME*)TAU,  (LAPACK_COMPLEX_TYPENAME*)&opt_work, &MONE, &info);
  #else
  SYTRD(uplo, N,  (LAPACK_COMPLEX_TYPENAME*)*aa, N, D, E,  (LAPACK_COMPLEX_TYPENAME*)TAU,  (LAPACK_COMPLEX_TYPENAME*)&opt_work, -1, &info);
  #endif
  if(info != 0) {
    if(a_out != NULL) free(a);
    return info;
  }
  int OWORK = (int)
      #ifdef LLW_real_double_IS_COMPLEX
              opt_work.re
      #else
              opt_work
            #endif
    ;
  LLW_real_double_COMPLEX* work = malloc(sizeof(LLW_real_double_COMPLEX)*(OWORK));
  #ifdef FCALL
  SYTRD(&uplo, &N,  (LAPACK_COMPLEX_TYPENAME*)*aa, &N, D, E,  (LAPACK_COMPLEX_TYPENAME*)TAU,  (LAPACK_COMPLEX_TYPENAME*)work, &OWORK, &info);
  #else
  SYTRD(uplo, N,  (LAPACK_COMPLEX_TYPENAME*)*aa, N, D, E,  (LAPACK_COMPLEX_TYPENAME*)TAU,  (LAPACK_COMPLEX_TYPENAME*)work, OWORK, &info);
  #endif
  #endif
  if(a_out != NULL) free(a);
  return info;
}

lapackwrap_real_double.c File Reference

Defines

Functions

Define Documentation

Function Documentation