lapackwrap_real_single.h File Reference

#include <liblapackwrap/lapackwrap_types_real_single.h>

Include dependency graph for lapackwrap_real_single.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Defines
#define	LLW_CONST   const
Functions
int	llw_real_single_stegr (char jobz, char range, int N, LLW_CONST LLW_real_single_REAL *D, LLW_CONST LLW_real_single_REAL *E, LLW_real_single_REAL VL, LLW_real_single_REAL VU, int IL, int IU, int *M, LLW_real_single_REAL *W, LLW_real_single_COMPLEX *Z, int LDZ, int *ISUPPZ)
int	llw_real_single_stegr_a (char jobz, char range, int N, LLW_CONST LLW_real_single_REAL *D, LLW_CONST LLW_real_single_REAL *E, LLW_real_single_REAL VL, LLW_real_single_REAL VU, int IL, int IU, int *M, LLW_real_single_REAL **W, LLW_real_single_COMPLEX **Z, int *LDZ, int **ISUPPZ)
int	llw_real_single_sbtrd_a (char vect, char uplo, int N, int KD, LLW_CONST LLW_real_single_COMPLEX *A, int LDAB, LLW_real_single_REAL **D, LLW_real_single_REAL **E, LLW_real_single_COMPLEX **Q)
int	llw_real_single_sbtrd (char vect, char uplo, int N, int KD, LLW_CONST LLW_real_single_COMPLEX *AB, int LDAB, LLW_real_single_REAL *D, LLW_real_single_REAL *E, LLW_real_single_COMPLEX *Q, int LDQ)
int	llw_real_single_sytrd (char uplo, int N, LLW_CONST LLW_real_single_COMPLEX *A, LLW_real_single_REAL *D, LLW_real_single_REAL *E, LLW_real_single_COMPLEX *TAU, LLW_real_single_COMPLEX **a_out)
int	llw_real_single_orgtr (char uplo, int N, LLW_real_single_COMPLEX *A, LLW_real_single_COMPLEX *TAU)

---

Define Documentation

#define LLW_CONST   const

Note:

Generic LLW data types are placed in this header file. defining LLW_real_single_FOO helps translate the generic core lapackwrap and blaswrap patterns to the final version. An alternative method would be translating each genericized core datatype into the specific data type. However, that would require a more complicated, multi-substitution method (e.g. with a perl script. This may come later, when I've gotten more experience using this lib.This is to ensure safety. If you don't care, you can save some cycles at the expense of "guaranteeing" safety for the data you're passing in read-only by commenting out this define. This will make all const arrays non-const and therefore not "guaranteed" to not be modified by the routines. This is effected, in addition to the const keyword, by copying each array into a secondary one via memcopy(). Turning this off saves that copy, which is the largst cost of this safety.

Definition at line 48 of file lapackwrap_real_single.h.

---

Function Documentation

int llw_real_single_orgtr	(	char	uplo,
		int	N,
		LLW_real_single_COMPLEX *	A,
		LLW_real_single_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_single.c.

References LAPACK_COMPLEX_TYPENAME, LLW_real_single_COMPLEX, and ORGTR.

                                                                                                      {
  int info;
  #if LAPACK == ACML
  ORGTR(uplo, N, (LAPACK_COMPLEX_TYPENAME*)A, N, (LAPACK_COMPLEX_TYPENAME*)TAU, &info);
  #else
  LLW_real_single_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_single_IS_COMPLEX
              opt_work.re
      #else
              opt_work
            #endif
    ;
  LLW_real_single_COMPLEX* work = (LLW_real_single_COMPLEX*)malloc(sizeof(LLW_real_single_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_single_sbtrd	(	char	vect,
		char	uplo,
		int	N,
		int	KD,
		LLW_CONST LLW_real_single_COMPLEX *	AB,
		int	LDAB,
		LLW_real_single_REAL *	D,
		LLW_real_single_REAL *	E,
		LLW_real_single_COMPLEX *	Q,
		int	LDQ
	)

Gets the tri-diagonal matrix for a banded matrix.

Definition at line 188 of file lapackwrap_real_single.c.

References LAPACK_COMPLEX_TYPENAME, LLW_real_single_COMPLEX, and SBTRD.

Referenced by llw_real_single_sbtrd_a().

                                                                                                                                                                                                       {
  LLW_real_single_COMPLEX *ab;
  #ifdef LLW_CONST
  ab = (LLW_real_single_COMPLEX*)malloc(N*(KD+1)*sizeof(LLW_real_single_COMPLEX));
  if(ab == NULL) return 100;
  memcpy(ab, AB, N*(KD+1)*sizeof(LLW_real_single_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_single_COMPLEX *work = (LLW_real_single_COMPLEX*)malloc(N*sizeof(LLW_real_single_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_single_sbtrd_a	(	char	vect,
		char	uplo,
		int	N,
		int	KD,
		LLW_CONST LLW_real_single_COMPLEX *	A,
		int	LDAB,
		LLW_real_single_REAL **	D,
		LLW_real_single_REAL **	E,
		LLW_real_single_COMPLEX **	Q
	)

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

Definition at line 171 of file lapackwrap_real_single.c.

References LLW_real_single_COMPLEX, LLW_real_single_REAL, and llw_real_single_sbtrd().

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

Here is the call graph for this function:

int llw_real_single_stegr	(	char	jobz,
		char	range,
		int	N,
		LLW_CONST LLW_real_single_REAL *	D,
		LLW_CONST LLW_real_single_REAL *	E,
		LLW_real_single_REAL	VL,
		LLW_real_single_REAL	VU,
		int	IL,
		int	IU,
		int *	M,
		LLW_real_single_REAL *	W,
		LLW_real_single_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_single.c.

References LAPACK_COMPLEX_TYPENAME, LLW_real_single_REAL, and STEGR.

Referenced by llw_real_single_stegr_a().

                                                                                                                                                                                                                                                                                   {
  /*
  fprintf(stderr, "llw_real_single_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_single_REAL *e_in;
  LLW_real_single_REAL *d_in;
  #ifdef LLW_CONST
  e_in = malloc(N*sizeof(LLW_real_single_REAL));
  d_in = malloc(N*sizeof(LLW_real_single_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_single_REAL));
  memcpy(d_in, D, N*sizeof(LLW_real_single_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_single_REAL opt_lwork;
  int opt_liwork;
#ifdef FCALL
  /*Look at all the fun in a FORTRAN call. Whee!*/
  LLW_real_single_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_single_REAL *work = malloc(sizeof(LLW_real_single_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_single_stegr_a	(	char	jobz,
		char	range,
		int	N,
		LLW_CONST LLW_real_single_REAL *	D,
		LLW_CONST LLW_real_single_REAL *	E,
		LLW_real_single_REAL	VL,
		LLW_real_single_REAL	VU,
		int	IL,
		int	IU,
		int *	M,
		LLW_real_single_REAL **	W,
		LLW_real_single_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_single_REAL; will be used to store eigvals
• Z pointer to a pointer to a LLW_real_single_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_single_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_single.c.

References LLW_real_single_COMPLEX, LLW_real_single_REAL, and llw_real_single_stegr().

                                                                                                                                                                                                                                                                                         {
  *W = malloc(sizeof(LLW_real_single_REAL)*N);
  if(*W == NULL) {
    return 100;
  }
  /*Right now, we're just doing N.*/
  *LDZ = N;
  *Z = malloc(sizeof(LLW_real_single_COMPLEX)*(*LDZ)*N);
  if(*Z == NULL) {
    free(*W);
    *W = NULL;
    return 100;
  }
  *ISUPPZ = malloc(sizeof(LLW_real_single_REAL)*2*N);
  if(*ISUPPZ == NULL) {
    free(*W);
    *W = NULL;
    free(*Z);
    *Z = NULL;
    return 100;
  }
  int info = llw_real_single_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_single_sytrd	(	char	uplo,
		int	N,
		LLW_CONST LLW_real_single_COMPLEX *	A,
		LLW_real_single_REAL *	D,
		LLW_real_single_REAL *	E,
		LLW_real_single_COMPLEX *	TAU,
		LLW_real_single_COMPLEX **	a_out
	)

Gets the tridiagonal matrix for a symmetric matrix.

Parameters:

uplo	'U' if upper triangle is desired, 'L' if lower triangle
N	size of an edge of A (since it's symmetric, is NxN)
A	original matrix
D	the diagonal of the final tridiag matrix
E	the upper (lower) subdiagonal.
TAU	LLW_CONSTs (needed for llw_real_single_orgr to get Q
a_out	pointer to a pointer to a LLW_real_single_REAL; used to save the resulting matrix (is A with output stored in it, as described by netlib's documentation on sytrd.

Returns:

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

Note:

A is LLW_CONST, so it will be preserved. If you want the tri-diagonal matrix, use a_out.

LAPACK wrapper for DSYTRD.

Returns:

the resulting INFO

Definition at line 223 of file lapackwrap_real_single.c.

References LAPACK_COMPLEX_TYPENAME, LLW_real_single_COMPLEX, and SYTRD.

                                                                                                                                                                                                   {
  LLW_real_single_COMPLEX **aa;
  LLW_real_single_COMPLEX *a;
  if(a_out != NULL) {
    a = *a_out;
  }else{
    /*Point aa to a LLW_real_single_COMPLEX pointer on the stack if the user doesn't want it*/
    aa = &a;
  }
  *aa = malloc(sizeof(LLW_real_single_COMPLEX)*N*N);
  if(a == NULL) return 100;
  memcpy(*aa, A, N*N*sizeof(LLW_real_single_COMPLEX));
  int info;
  #if LAPACK == ACML
  SYTRD(uplo, N,  (LAPACK_COMPLEX_TYPENAME*)*aa, N, D, E,  (LAPACK_COMPLEX_TYPENAME*)TAU, &info);
  #else
  LLW_real_single_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_single_IS_COMPLEX
              opt_work.re
      #else
              opt_work
            #endif
    ;
  LLW_real_single_COMPLEX* work = malloc(sizeof(LLW_real_single_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;
}