liblapackwrap 0.1
Defines | Functions

lapackwrap_complex_double.h File Reference

#include <liblapackwrap/lapackwrap_types_complex_double.h>
Include dependency graph for lapackwrap_complex_double.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_complex_double_stegr (char jobz, char range, int N, LLW_CONST LLW_complex_double_REAL *D, LLW_CONST LLW_complex_double_REAL *E, LLW_complex_double_REAL VL, LLW_complex_double_REAL VU, int IL, int IU, int *M, LLW_complex_double_REAL *W, LLW_complex_double_COMPLEX *Z, int LDZ, int *ISUPPZ)
int llw_complex_double_stegr_a (char jobz, char range, int N, LLW_CONST LLW_complex_double_REAL *D, LLW_CONST LLW_complex_double_REAL *E, LLW_complex_double_REAL VL, LLW_complex_double_REAL VU, int IL, int IU, int *M, LLW_complex_double_REAL **W, LLW_complex_double_COMPLEX **Z, int *LDZ, int **ISUPPZ)
int llw_complex_double_sbtrd_a (char vect, char uplo, int N, int KD, LLW_CONST LLW_complex_double_COMPLEX *A, int LDAB, LLW_complex_double_REAL **D, LLW_complex_double_REAL **E, LLW_complex_double_COMPLEX **Q)
int llw_complex_double_sbtrd (char vect, char uplo, int N, int KD, LLW_CONST LLW_complex_double_COMPLEX *AB, int LDAB, LLW_complex_double_REAL *D, LLW_complex_double_REAL *E, LLW_complex_double_COMPLEX *Q, int LDQ)
int llw_complex_double_sytrd (char uplo, int N, LLW_CONST LLW_complex_double_COMPLEX *A, LLW_complex_double_REAL *D, LLW_complex_double_REAL *E, LLW_complex_double_COMPLEX *TAU, LLW_complex_double_COMPLEX **a_out)
int llw_complex_double_orgtr (char uplo, int N, LLW_complex_double_COMPLEX *A, LLW_complex_double_COMPLEX *TAU)

Define Documentation

#define LLW_CONST   const
Note:
Generic LLW data types are placed in this header file. defining LLW_complex_double_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_complex_double.h.

Function Documentation

int llw_complex_double_orgtr ( char  uplo,
int  N,
LLW_complex_double_COMPLEX *  A,
LLW_complex_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_complex_double.c.

References LAPACK_COMPLEX_TYPENAME, LLW_complex_double_COMPLEX, LLW_complex_double_IS_COMPLEX, and ORGTR.

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

Gets the tri-diagonal matrix for a banded matrix.

Definition at line 188 of file lapackwrap_complex_double.c.

References LAPACK_COMPLEX_TYPENAME, LLW_complex_double_COMPLEX, and SBTRD.

Referenced by llw_complex_double_sbtrd_a().

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

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

Definition at line 171 of file lapackwrap_complex_double.c.

References LLW_complex_double_COMPLEX, LLW_complex_double_REAL, and llw_complex_double_sbtrd().

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

Here is the call graph for this function:

int llw_complex_double_stegr ( char  jobz,
char  range,
int  N,
LLW_CONST LLW_complex_double_REAL *  D,
LLW_CONST LLW_complex_double_REAL *  E,
LLW_complex_double_REAL  VL,
LLW_complex_double_REAL  VU,
int  IL,
int  IU,
int *  M,
LLW_complex_double_REAL *  W,
LLW_complex_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_complex_double.c.

References LAPACK_COMPLEX_TYPENAME, LLW_complex_double_REAL, and STEGR.

Referenced by llw_complex_double_stegr_a().

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

References LLW_complex_double_COMPLEX, LLW_complex_double_REAL, and llw_complex_double_stegr().

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

Gets the tridiagonal matrix for a symmetric matrix.

Parameters:
uplo'U' if upper triangle is desired, 'L' if lower triangle
Nsize of an edge of A (since it's symmetric, is NxN)
Aoriginal matrix
Dthe diagonal of the final tridiag matrix
Ethe upper (lower) subdiagonal.
TAULLW_CONSTs (needed for llw_complex_double_orgr to get Q
a_outpointer to a pointer to a LLW_complex_double_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_complex_double.c.

References LAPACK_COMPLEX_TYPENAME, LLW_complex_double_COMPLEX, LLW_complex_double_IS_COMPLEX, and SYTRD.

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