evsl_f90.c

Go to the documentation of this file.

#include <string.h>
#include <stdint.h>
#include "def.h"
#include "struct.h"
#include "internal_proto.h"

int evsl_nev_computed=0, evsl_n=0;
double *evsl_eigval_computed=NULL, *evsl_eigvec_computed=NULL;

void EVSLFORT(evsl_start)() {
  EVSLStart();
}

void EVSLFORT(evsl_finish)() {
  EVSLFinish();
}

void EVSLFORT(evsl_coo2csr)(int *n, int *nnz, int *ir, int *jc, 
                            double *vv, uintptr_t *csrf90) {
  cooMat coo;
  coo.nrows = *n;
  coo.ncols = *n;
  coo.nnz = *nnz;
  coo.ir = ir;
  coo.jc = jc;
  coo.vv = vv;
  csrMat *csr;
  Malloc(csr, 1, csrMat);
  cooMat_to_csrMat(1, &coo, csr);
  *csrf90 = (uintptr_t) csr;
}

void EVSLFORT(evsl_arr2csr)(int *n, int *ia, int *ja, 
                            double *a, uintptr_t *csrf90) {
  csrMat *csr;
  Malloc(csr, 1, csrMat);
  /* does not own the data */
  csr->owndata = 0;
  csr->nrows = *n;
  csr->ncols = *n;
  csr->ia = ia;
  csr->ja = ja;
  csr->a = a;
  *csrf90 = (uintptr_t) csr;
}

void EVSLFORT(evsl_free_csr)(uintptr_t *csrf90) {
  csrMat *csr = (csrMat *) (*csrf90);
  free_csr(csr);
  free(csr);
}

void EVSLFORT(evsl_seta_csr)(uintptr_t *Af90) {
  csrMat *A = (csrMat *) (*Af90);
  SetAMatrix(A);
}

void EVSLFORT(evsl_setb_csr)(uintptr_t *Bf90) {
  csrMat *B = (csrMat *) (*Bf90);
  SetBMatrix(B);
}

void EVSLFORT(evsl_setamv)(int *n, void *func, void *data) {
  SetAMatvec(*n, (MVFunc) func, data);
}

void EVSLFORT(evsl_setbmv)(int *n, void *func, void *data) {
  SetBMatvec(*n, (MVFunc) func, data);
}

void EVSLFORT(evsl_setbsol)(void *func, void *data) {
  SetBSol((SolFuncR) func, data);
}

void EVSLFORT(evsl_setltsol)(void *func, void *data) {
  SetLTSol((SolFuncR) func, data);
}

void EVSLFORT(evsl_set_geneig)() {
  SetGenEig();
}

void EVSLFORT(evsl_lanbounds)(int *nsteps, double *lmin, double *lmax) {
  int n = evsldata.n;
  double *vinit;
  Malloc(vinit, n, double);
  rand_double(n, vinit);
  LanTrbounds(50, *nsteps, 1e-10, vinit, 1, lmin, lmax, NULL);
  //LanBounds(*nsteps, vinit, lmin, lmax);
  free(vinit);
}

void EVSLFORT(evsl_kpm_spslicer)(int *Mdeg, int *nvec, double *xintv,
                                int *nslices, double *sli, int *evint) {
  int npts;
  double *mu, ecount;
  Malloc(mu, *Mdeg+1, double);
  kpmdos(*Mdeg, 1, *nvec, xintv, mu, &ecount);
  fprintf(stdout, " estimated eig count in interval: %.15e \n",ecount);
  
  npts = 10 *ecount;
  int ierr = spslicer(sli, mu, *Mdeg, xintv, *nslices, npts);
  if (ierr) {
    printf("spslicer error %d\n", ierr);
    exit(1);
  }
  *evint = (int) (1 + ecount / (*nslices));
  free(mu);
}

void EVSLFORT(evsl_find_pol)(double *xintv, double *thresh_int, 
                             double *thresh_ext, uintptr_t *polf90) {
  polparams *pol;
  Malloc(pol, 1, polparams);
  set_pol_def(pol);
  pol->damping = 2;
  pol->thresh_int = *thresh_int;
  pol->thresh_ext = *thresh_ext;
  pol->max_deg  = 300;
  find_pol(xintv, pol);       
  
  fprintf(stdout, " polynomial deg %d, bar %e gam %e\n",
          pol->deg,pol->bar, pol->gam);

  *polf90 = (uintptr_t) pol;
}

void EVSLFORT(evsl_free_pol)(uintptr_t *polf90) {
  /* cast pointer */
  polparams *pol = (polparams *) (*polf90);
  free_pol(pol);
  free(pol);
}

void EVSLFORT(evsl_find_rat)(double *intv, uintptr_t *ratf90) {
  int pow = 2;
  int num = 1;
  double beta = 0.01;
  ratparams *rat;
  Malloc(rat, 1, ratparams);
  set_ratf_def(rat);
  rat->pw = pow;
  rat->num = num;
  rat->beta = beta;
  find_ratf(intv, rat);

  *ratf90 = (uintptr_t) rat;
}

void EVSLFORT(evsl_free_rat)(uintptr_t *ratf90) {
  ratparams *rat = (ratparams *) (*ratf90);
  free_rat(rat);
  free(rat);
}

void EVSLFORT(evsl_cheblantr)(int *mlan, int *nev, double *xintv, int *max_its,
                              double *tol, uintptr_t *polf90) {
  int n, nev2, ierr;
  double *lam, *Y, *res;
  FILE *fstats = stdout;
  double *vinit;
 
  n = evsldata.n;
  Malloc(vinit, n, double);
  rand_double(n, vinit);

  /* cast pointer */
  polparams *pol = (polparams *) (*polf90);

  ierr = ChebLanTr(*mlan, *nev, xintv, *max_its, *tol, vinit,
                   pol, &nev2, &lam, &Y, &res, fstats);

  if (ierr) {
    printf("ChebLanTr error %d\n", ierr);
  }

  free(vinit);
  if (res) {
    free(res);
  }
  /* save pointers to the global variables */
  evsl_nev_computed = nev2;
  evsl_n = n;
  evsl_eigval_computed = lam;
  evsl_eigvec_computed = Y;
}

void EVSLFORT(evsl_cheblannr)(double *xintv, int *max_its, double *tol, 
                              uintptr_t *polf90) {
  int n, nev2, ierr;
  double *lam, *Y, *res;
  FILE *fstats = stdout;
  double *vinit;
 
  n = evsldata.n;
  Malloc(vinit, n, double);
  rand_double(n, vinit);

  /* cast pointer */
  polparams *pol = (polparams *) (*polf90);

  ierr = ChebLanNr(xintv, *max_its, *tol, vinit,
                   pol, &nev2, &lam, &Y, &res, fstats);

  if (ierr) {
    printf("ChebLanNr error %d\n", ierr);
  }

  free(vinit);
  if (res) {
    free(res);
  }
  /* save pointers to the global variables */
  evsl_nev_computed = nev2;
  evsl_n = n;
  evsl_eigval_computed = lam;
  evsl_eigvec_computed = Y;
}

void EVSLFORT(evsl_ratlannr)(double *xintv, int *max_its, double *tol, 
                             uintptr_t *ratf90) {
  int n, nev2, ierr;
  double *lam, *Y, *res;
  FILE *fstats = stdout;
  double *vinit;
 
  n = evsldata.n;
  Malloc(vinit, n, double);
  rand_double(n, vinit);

  /* cast pointer */
  ratparams *rat = (ratparams *) (*ratf90);

  ierr = RatLanNr(xintv, *max_its, *tol, vinit,
                  rat, &nev2, &lam, &Y, &res, fstats);

  if (ierr) {
    printf("RatLanNr error %d\n", ierr);
  }

  free(vinit);
  if (res) {
    free(res);
  }
  /* save pointers to the global variables */
  evsl_nev_computed = nev2;
  evsl_n = n;
  evsl_eigval_computed = lam;
  evsl_eigvec_computed = Y;
}

void EVSLFORT(evsl_ratlantr)(int * lanm, int * nev, double *xintv, 
        int *max_its, double *tol, uintptr_t *ratf90) {
  int n, nev2, ierr;
  double *lam, *Y, *res;
  FILE *fstats = stdout;
  double *vinit;
 
  n = evsldata.n;
  Malloc(vinit, n, double);
  rand_double(n, vinit);

  /* cast pointer */
  ratparams *rat = (ratparams *) (*ratf90);

  ierr = RatLanTr(*lanm, *nev, xintv, *max_its, *tol, vinit,
                  rat, &nev2, &lam, &Y, &res, fstats);

  if (ierr) {
    printf("RatLanNr error %d\n", ierr);
  }

  free(vinit);
  if (res) {
    free(res);
  }
  /* save pointers to the global variables */
  evsl_nev_computed = nev2;
  evsl_n = n;
  evsl_eigval_computed = lam;
  evsl_eigvec_computed = Y;
}

void EVSLFORT(evsl_get_nev)(int *nev) {
  *nev = evsl_nev_computed;
}

void EVSLFORT(evsl_copy_result)(double *val, double *vec) {
  memcpy(val, evsl_eigval_computed, evsl_nev_computed*sizeof(double));
  memcpy(vec, evsl_eigvec_computed, evsl_nev_computed*evsl_n*sizeof(double));
  /* reset global variables */
  evsl_nev_computed = 0;
  evsl_n = 0;
  free(evsl_eigval_computed);
  free(evsl_eigvec_computed);
  evsl_eigval_computed = NULL;
  evsl_eigvec_computed = NULL;
}