lapackwrap_real_double.c

Go to the documentation of this file.

/*
 * Wrapper for a variety of LAPACK routines, to provide a unified interface.

 Copyright (C) 2010 Joseph Pingenot

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

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

#define ACML 0
#define CELL 1

/*These follow the FORTRAN calling semantics.*/
#if LAPACK == CELL
#define FCALL
#endif

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) {
  *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;
}

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) {
  /*
  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;
}



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) {
  *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);
}

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) {
  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;
}



/**LAPACK wrapper for DSYTRD.
 *\return the resulting INFO
 */
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) {
  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;
}



int llw_real_double_orgtr(char uplo, int N, LLW_real_double_COMPLEX* A, LLW_real_double_COMPLEX* TAU) {
  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;
}