evsl_suitesparse.c

Go to the documentation of this file.

#include <string.h>
#include "def.h"
#include "struct.h"
#include "internal_proto.h"
#include "cholmod.h"
#include "umfpack.h"
#include "evsl_direct.h"

#ifdef CHOLMOD_USE_LONG
#define CHOLMOD(name) cholmod_l_ ## name
#else
#define CHOLMOD(name) cholmod_ ## name
#endif

typedef struct _BSolDataDirect {
  cholmod_factor *LB;
  cholmod_common cm;
} BSolDataDirect;

cholmod_sparse* csrMat_to_cholmod_sparse(csrMat *A, int stype) {
  cholmod_sparse *B = NULL;
  Malloc(B, 1, cholmod_sparse);
  B->nrow = A->nrows;
  B->ncol = A->ncols;
  B->nzmax = A->ia[A->nrows];
#ifdef CHOLMOD_USE_LONG
  /* long int version, allocate memory and copy ia, ja */
  Malloc(B->p, B->nrow+1, SuiteSparse_long);
  Malloc(B->i, B->nzmax, SuiteSparse_long);
  int i;
  for (i=0; i<B->nrow+1; i++) {
    SuiteSparse_long *ptr = (SuiteSparse_long *) B->p;
    ptr[i] = A->ia[i];
  }
  for (i=0; i<B->nzmax; i++) {
    SuiteSparse_long *ptr = (SuiteSparse_long *) B->i;
    ptr[i] = A->ja[i];
  }
#else
  /* column pointers */
  B->p = A->ia;
  /* row indices */
  B->i = A->ja;
#endif
  B->nz = NULL;
  B->x = A->a;
  B->z = NULL;
  B->stype = stype;
#ifdef CHOLMOD_USE_LONG
  B->itype = CHOLMOD_LONG;
#else
  B->itype = CHOLMOD_INT;
#endif
  B->xtype = CHOLMOD_REAL;
  B->dtype = CHOLMOD_DOUBLE;
  B->sorted = 0;
  B->packed = 1;

  return B;
}

void arr_to_cholmod_dense(int m, int n, double *A, cholmod_dense *B) {
  B->nrow = m;
  B->ncol = n;
  B->nzmax = m*n;
  B->d = m;
  B->x = A;
  B->z = NULL;
  B->xtype = CHOLMOD_REAL;
  B->dtype = CHOLMOD_DOUBLE;
}

#if 0

cholmod_dense* arr_to_cholmod_dense_alloc(int m, int n, double *A) {
  cholmod_dense *B = NULL;
  Malloc(B, 1, cholmod_dense);
  B->nrow = m;
  B->ncol = n;
  B->nzmax = m*n;
  B->d = m;
  B->x = A;
  B->z = NULL;
  B->xtype = CHOLMOD_REAL;
  B->dtype = CHOLMOD_DOUBLE;

  return B;
}
#endif


cholmod_dense cholmod_X, cholmod_B, *cholmod_Y=NULL, *cholmod_E=NULL,
              *cholmod_W=NULL;

void BSolDirect(double *b, double *x, void *data) {
  //int n;

  BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
  cholmod_factor *LB;
  cholmod_common *cc;

  LB = Bsol_data->LB;
  cc = &Bsol_data->cm;
  //n = LB->n;

  /* give the wrapper data */
  cholmod_X.x = x;
  cholmod_B.x = b;

  /* XXX is this always safe for cholmod_x ? */
  cholmod_dense *cholmod_Xp = &cholmod_X; 
  CHOLMOD(solve2)(CHOLMOD_A, LB, &cholmod_B, NULL, 
                  &cholmod_Xp, NULL, &cholmod_Y, &cholmod_E, cc);
}

int SetupBSolDirect(csrMat *B, void **data) {
  double tms = evsl_timer();
  BSolDataDirect *Bsol_data;
  Malloc(Bsol_data, 1, BSolDataDirect);
  cholmod_sparse *Bcholmod;
  cholmod_common *cc = &Bsol_data->cm;
  cholmod_factor *LB;

  /* start CHOLMOD */
  CHOLMOD(start)(cc);
  /* force to have LL factor */
  cc->final_asis = 0;
  cc->final_ll = 1;
  /* convert matrix. 
   * stype=1 means the upper triangular part of B will be accessed */
  Bcholmod = csrMat_to_cholmod_sparse(B, 1);
  /* check common and the matrix */
  CHOLMOD(check_common)(cc);
  CHOLMOD(check_sparse)(Bcholmod, cc);
  /* symbolic and numeric fact */
  //printf("start analysis\n");
  LB = CHOLMOD(analyze)(Bcholmod, cc);
  CHOLMOD(check_common)(cc);
  //printf("cc status %d\n", cc->status);
  //printf("start factorization\n");
  CHOLMOD(factorize)(Bcholmod, LB, cc);
  //printf("done factorization\n");
  /* check the factor */
  CHOLMOD(check_factor)(LB, cc);
  /* save the struct to global variable */
  Bsol_data->LB = LB;
  /* check the factor */
  CHKERR(LB->is_ll == 0);
  /* free the matrix data and wrapper */
#ifdef CHOLMOD_USE_LONG
  free(Bcholmod->p);
  free(Bcholmod->i);
#endif
  /* free the wrapper itself */
  free(Bcholmod);
  /* space for solve */
  int n = B->nrows;
  /* setup the wrapper, no data given yet */
  arr_to_cholmod_dense(n, 1, NULL, &cholmod_X);
  arr_to_cholmod_dense(n, 1, NULL, &cholmod_B);

  *data = (void *) Bsol_data;

  double tme = evsl_timer();
  evslstat.t_setBsv += tme - tms;

  return 0;
}

void FreeBSolDirectData(void *vdata) {
  BSolDataDirect *data = (BSolDataDirect *) vdata;
  cholmod_common *cc = &data->cm;
  /* free the factor */
  CHOLMOD(free_factor)(&data->LB, cc);
  /* free workspace for solve2 */ 
  CHOLMOD(free_dense)(&cholmod_Y, cc);
  CHOLMOD(free_dense)(&cholmod_E, cc);
  CHOLMOD(free_dense)(&cholmod_W, cc);
  /* finish cholmod */
  CHOLMOD(finish)(cc);
  free(vdata);
}

void LTSolDirect(double *b, double *x, void *data) {
  //int n;

  BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
  cholmod_factor *LB;
  cholmod_common *cc;

  LB = Bsol_data->LB;
  cc = &Bsol_data->cm;
  //n = LB->n;
  
  /* XXX are these always safe ? */
  cholmod_B.x = b;
  CHOLMOD(solve2)(CHOLMOD_Lt, LB, &cholmod_B, NULL, 
                  &cholmod_W, NULL, &cholmod_Y, &cholmod_E, cc);
  
  cholmod_X.x = x;
  cholmod_dense *cholmod_Xp = &cholmod_X;
  CHOLMOD(solve2)(CHOLMOD_Pt, LB, cholmod_W, NULL, 
                  &cholmod_Xp, NULL, &cholmod_Y, &cholmod_E, cc);
}

void ASIGMABSolDirect(int n, double *br, double *bi, double *xr, 
                      double *xz, void *data) {
  void* Numeric = data;
  double Control[UMFPACK_CONTROL]; 
  umfpack_zl_defaults(Control);
  Control[UMFPACK_IRSTEP] = 0; // no iterative refinement for umfpack
  umfpack_zl_solve(UMFPACK_A, NULL, NULL, NULL, NULL, xr, xz, br, bi, 
                   Numeric, Control, NULL);
}

int SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num,
                          complex double *zk, void **data) {
  int i, j, nrow, ncol, nnzB, nnzC, *map, status;
  csrMat *B, C, eye;
  double tms = evsl_timer();
  /* UMFPACK matrix for the shifted matrix 
   * C = A - s * B */
  SuiteSparse_long *Cp, *Ci;
  double *Cx, *Cz, zkr1;
  void *Symbolic=NULL, *Numeric=NULL;
  
  nrow = A->nrows;
  ncol = A->ncols;

  if (BB) {
    B = BB;
  } else {
    /* if B==NULL, B=I, standard e.v. prob */
    speye(nrow, &eye);
    B = &eye;
  }

  nnzB = B->ia[nrow];
  /* map from nnz in B to nnz in C, useful for multi-poles */
  Malloc(map, nnzB, int);
  /* NOTE: SuiteSparse requires that the matrix entries be sorted.
   * The matadd routine will  guarantee this
   * C = A + 0.0 * B 
   * This actually computes the pattern of A + B
   * and also can guarantee C has sorted rows 
   * map is the mapping from nnzB to nnzC */
  matadd(1.0, 0.0, A, B, &C, NULL, map);
  nnzC = C.ia[nrow];
  /* malloc and copy to SuiteSparse matrix */
  Malloc(Cp, nrow+1, SuiteSparse_long);
  Malloc(Ci, nnzC, SuiteSparse_long);
  Calloc(Cz, nnzC, double);
  for (i=0; i<nrow+1; i++) {
    Cp[i] = C.ia[i];
  }
  for (i=0; i<nnzC; i++) {
    Ci[i] = C.ja[i];
  }
  Cx = C.a;
  /* pole loop
   * for each pole we shift with B and factorize */
  zkr1 = 0.0;
  for (i=0; i<num; i++) {
    /* the complex shift for pole i */
    double zkr = creal(zk[i]);
    double zkc = cimag(zk[i]);

    // shift B
    for (j=0; j<nnzB; j++) {
      int p = map[j];
      double v = B->a[j];
      CHKERR(Ci[p] != B->ja[j]);
      Cx[p] -= (zkr - zkr1) * v;
      Cz[p] = -zkc * v;
    }
    /* only do symbolic factorization once */
    if (i==0) {
      /* Symbolic Factorization */
      status = umfpack_zl_symbolic(nrow, ncol, Cp, Ci, Cx, Cz, &Symbolic, 
                                   NULL, NULL);
      if (status < 0) {
        printf("umfpack_zl_symbolic failed, %d\n", status);
        return 1;
      }
    }
    /* Numerical Factorization */
    double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO];
    Control[UMFPACK_PRL] = 0;
    status = umfpack_zl_numeric(Cp, Ci, Cx, Cz, Symbolic, &Numeric, Control, Info);
    umfpack_zl_report_info(Control, Info);
    if (status < 0) {
      printf("umfpack_zl_numeric failed and exit, %d\n", status);
      return 1;
    }
    /* save the data */
    data[i] = Numeric;
    /* for the next shift */
    zkr1 = zkr;
  } /* for (i=...)*/
  
  /* free the symbolic fact */
  if (Symbolic) {
    umfpack_zl_free_symbolic(&Symbolic);
  }

  free(map);
  free(Cp);
  free(Ci);
  free(Cz);
  free_csr(&C);
  if (!BB) {
    free_csr(&eye);
  }

  double tme = evsl_timer();
  evslstat.t_setASigBsv += tme - tms;
  
  return 0;
}

void FreeASIGMABSolDirect(int num, void **data) {
  int i;
  for (i=0; i<num; i++) {
    umfpack_zl_free_numeric(&data[i]);
  }
}