blaswrap_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/blaswrap_real_double.h"
#include "liblapackwrap/blaswrap_functions_real_double.c"

/**frontend for llw_real_double_gemm but if you don't have a matrix C yet (clearly, if you're calling this routine, you want one.
 */
int llw_real_double_gemm_a(char transa, char transb, int M, int N, int K, LLW_real_double_COMPLEX *alpha, LLW_CONST LLW_real_double_COMPLEX *A, int LDA, LLW_CONST LLW_real_double_COMPLEX *B, int LDB, LLW_real_double_COMPLEX **C) {
  *C = (LLW_real_double_COMPLEX*)malloc(M*N*sizeof(LLW_real_double_COMPLEX));
  if((*C) == NULL) return 100;
  //fprintf(stderr, "N=%d M=%d N*M=%lu\n", M, N, (long unsigned)M*N);
  for(long unsigned i=0; i<M*N; i++) {
    //fprintf(stderr, "\tC=%p *C=%p (*C)[%lu]=%g\n", (void*)C, (void*)*C, i, (*C)[i]);
    llw_bw_real_double_set_from_double(&((*C)[i]), 0);
    //fprintf(stderr, "\t\t->(*C)[%lu]=%g\n", i, (*C)[i]);
  }
  LLW_real_double_COMPLEX zero;
  llw_bw_real_double_set_from_double(&zero, 0.0);
  return llw_real_double_gemm(transa, transb, M, N, K, alpha, A, LDA, B, LDB, &zero, *C, M);
}

int llw_real_double_gemm(char transa, char transb, int M, int N, int K, LLW_real_double_COMPLEX *alpha, LLW_CONST LLW_real_double_COMPLEX *A, int LDA, LLW_CONST LLW_real_double_COMPLEX *B, int LDB, LLW_real_double_COMPLEX *beta, LLW_real_double_COMPLEX *C, int LDC) {
  int ka=(transa=='N')?K:M;
  int kb=(transb=='N')?K:N;
  LLW_real_double_COMPLEX *a;
  LLW_real_double_COMPLEX *b;
  #ifdef LLW_CONST
  a = (LLW_real_double_COMPLEX*)malloc(LDA*ka*sizeof(LLW_real_double_COMPLEX));
  if(a == NULL) return 100;
  b = (LLW_real_double_COMPLEX*)malloc(LDA*kb*sizeof(LLW_real_double_COMPLEX));
  if(b == NULL) {
    free(a);
    return 100;
  }
  memcpy(a, A, LDA*ka*sizeof(LLW_real_double_COMPLEX));
  memcpy(b, B, LDB*kb*sizeof(LLW_real_double_COMPLEX));
  #else
  a = A;
  b = B;
  #endif
  #ifdef FCALL
  GEMM(&transa, &transb, &M, &N, &K, 
       #ifdef LLW_real_double_IS_COMPLEX
       (LAPACK_COMPLEX_TYPENAME*)alpha, 
       #else
       (LAPACK_COMPLEX_TYPENAME*)&alpha, 
       #endif
       (LAPACK_COMPLEX_TYPENAME*)a, &LDA, (LAPACK_COMPLEX_TYPENAME*)b, &LDB, 
       #ifdef LLW_real_double_IS_COMPLEX
       (LAPACK_COMPLEX_TYPENAME*)beta, 
       #else
       (LAPACK_COMPLEX_TYPENAME*)&beta, 
       #endif
       (LAPACK_COMPLEX_TYPENAME*)C, &LDC);
  #else
  GEMM(transa, transb, M, N, K, 
       #ifdef LLW_real_double_IS_COMPLEX
       (LAPACK_COMPLEX_TYPENAME*)alpha,
       #else
       *alpha,
       #endif
       (LAPACK_COMPLEX_TYPENAME*)a, LDA, (LAPACK_COMPLEX_TYPENAME*)b, LDB,
       #ifdef LLW_real_double_IS_COMPLEX
       (LAPACK_COMPLEX_TYPENAME*)beta,
       #else
       *beta,
       #endif
       (LAPACK_COMPLEX_TYPENAME*)C, LDC);
  #endif
  #ifdef LLW_CONST
  free(a);
  free(b);
  #endif
  return 0;
}