EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
lapl.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <float.h>
3 #include "def.h"
4 #include "struct.h"
5 #include "internal_proto.h"
6 
17 int lapgen(int nx, int ny, int nz, cooMat *Acoo) {
18  int n = nx * ny * nz;
19  Acoo->nrows = n;
20  Acoo->ncols = n;
21 
22  int nzmax = nz > 1 ? 7*n : 5*n;
23  Malloc(Acoo->ir, nzmax, int);
24  Malloc(Acoo->jc, nzmax, int);
25  Malloc(Acoo->vv, nzmax, double);
26 
27  int ii, nnz=0;
28  for (ii=0; ii<n; ii++) {
29  double v = -1.0;
30  int i,j,k,jj;
31  k = ii / (nx*ny);
32  i = (ii - k*nx*ny) / nx;
33  j = ii - k*nx*ny - i*nx;
34 
35  if (k > 0) {
36  jj = ii - nx * ny;
37  Acoo->ir[nnz] = ii; Acoo->jc[nnz] = jj; Acoo->vv[nnz] = v; nnz++;
38  }
39  if (k < nz-1) {
40  jj = ii + nx * ny;
41  Acoo->ir[nnz] = ii; Acoo->jc[nnz] = jj; Acoo->vv[nnz] = v; nnz++;
42  }
43 
44  if (i > 0) {
45  jj = ii - nx;
46  Acoo->ir[nnz] = ii; Acoo->jc[nnz] = jj; Acoo->vv[nnz] = v; nnz++;
47  }
48  if (i < ny-1) {
49  jj = ii + nx;
50  Acoo->ir[nnz] = ii; Acoo->jc[nnz] = jj; Acoo->vv[nnz] = v; nnz++;
51  }
52 
53  if (j > 0) {
54  jj = ii - 1;
55  Acoo->ir[nnz] = ii; Acoo->jc[nnz] = jj; Acoo->vv[nnz] = v; nnz++;
56  }
57  if (j < nx-1) {
58  jj = ii + 1;
59  Acoo->ir[nnz] = ii; Acoo->jc[nnz] = jj; Acoo->vv[nnz] = v; nnz++;
60  }
61 
62  v = nz > 1 ? 6.0 : 4.0;
63  Acoo->ir[nnz] = ii; Acoo->jc[nnz] = ii; Acoo->vv[nnz] = v; nnz++;
64  }
65 
66  Acoo->nnz = nnz;
67 
68  return 0;
69 }
70 
83 int exeiglap3(int nx, int ny, int nz, double a, double b, int *m, double **vo) {
84  double thetax = 0.5 * PI / (nx + 1.0);
85  double thetay = 0.5 * PI / (ny + 1.0);
86  double thetaz = 0.5 * PI / (nz + 1.0);
87  int i, j, k, l=0, n=nx*ny*nz;
88  double *v;
89  double tx, ty, tz, ev;
90  Malloc(v, n, double);
91  for (i=1; i<=nx; i++) {
92  tx = sin(i*thetax);
93  for (j=1; j<=ny; j++) {
94  ty = sin(j*thetay);
95  for (k=1; k<=nz; k++) {
96  tz = sin(k*thetaz);
97  if (1 == nz) {
98  tz = 0.0;
99  }
100  ev = 4*(tx*tx+ty*ty+tz*tz);
101  if (ev >= a - DBL_EPSILON && ev <= b + DBL_EPSILON) {
102  v[l++] = ev;
103  }
104  }
105  }
106  }
107  Realloc(v, l, double);
108  sort_double(l, v, NULL);
109 
110  *m = l;
111  *vo = v;
112  return 0;
113 }
114 
defs in EVSL
int nnz
Definition: struct.h:17
double * vv
Definition: struct.h:22
int * jc
Definition: struct.h:17
int * ir
Definition: struct.h:17
void sort_double(int n, double *v, int *ind)
Definition: vect.c:92
int exeiglap3(int nx, int ny, int nz, double a, double b, int *m, double **vo)
Exact eigenvalues of Laplacean in interval [a b].
Definition: lapl.c:83
This file contains function prototypes and constant definitions internally used in EVSL...
#define Malloc(base, nmem, type)
Definition: def.h:22
structs used in evsl
int ncols
Definition: struct.h:17
int lapgen(int nx, int ny, int nz, cooMat *Acoo)
Laplacean Matrix generator.
Definition: lapl.c:17
#define PI
Definition: def.h:15
int nrows
Definition: struct.h:17
#define Realloc(base, nmem, type)
Definition: def.h:42
sparse matrix format: the coordinate (COO) format, 0-based
Definition: struct.h:16