Skip to content

Commit

Permalink
Multi-dimensional differentiation working.
Browse files Browse the repository at this point in the history 
darcs-hash:20080228172636-40437-03d4a84297dec087eb9703dbf059cd967703cb07.gz
  • Loading branch information
@jedbrown
jedbrown committed Feb 28, 2008
1 parent ca83293 commit b2d95fb
Show file tree
Hide file tree
Showing 3 changed files with 211 additions and 14 deletions.
57 changes: 46 additions & 11 deletions cheb.c
View file
Expand Up @@ -17,39 +17,74 @@ int main(int argc,char **args)
{
Vec x,b,u; /* approx solution, RHS, exact solution */
Mat A;
Vec x2,b2,u2; /* approx solution, RHS, exact solution */
Mat A2;
KSP ksp; /* linear solver context */
PC pc; /* preconditioner context */
PetscReal norm; /* norm of solution error */
PetscScalar v,one = 1.0,none = -1.0;
PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,its, mn;
PetscInt i,j,Ii,J,Istart,Iend, m1 = 5, m = 8,n = 7, p=1, d=0, its;
PetscErrorCode ierr;
PetscTruth user_defined_pc, user_defined_mat;

ierr = PetscInitialize(&argc,&args,(char *)0,help); CHKERRQ(ierr);

ierr = fftw_import_system_wisdom(); CHKERRQ(ierr);

ierr = PetscOptionsGetInt(PETSC_NULL,"-m1",&m1,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);

ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,m);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-d",&d,PETSC_NULL);CHKERRQ(ierr);

// ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
// ierr = VecSetSizes(u,PETSC_DECIDE,m);CHKERRQ(ierr);
// ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = VecCreateSeq(PETSC_COMM_WORLD, m1, &u); CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);

ierr = MatCreateChebD1(PETSC_COMM_WORLD, x, b,
FFTW_ESTIMATE | FFTW_PRESERVE_INPUT, &A);CHKERRQ(ierr);
FFTW_ESTIMATE | FFTW_PRESERVE_INPUT, &A); CHKERRQ(ierr);

int dims[] = { m, n, p };
// ierr = VecCreate(PETSC_COMM_WORLD, &u2);CHKERRQ(ierr);
// ierr = VecSetSizes(u2, PETSC_DECIDE, m * n);CHKERRQ(ierr);
// ierr = VecSetFromOptions(u2);CHKERRQ(ierr);
ierr = VecCreateSeq(PETSC_COMM_WORLD, m * n * p, &u2); CHKERRQ(ierr);
ierr = VecDuplicate(u2, &b2);CHKERRQ(ierr);
ierr = VecDuplicate(b2, &x2);CHKERRQ(ierr);
ierr = MatCreateCheb(PETSC_COMM_WORLD, 3, d, dims, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT,
x2, b2, &A2); CHKERRQ(ierr);

ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

// Solution function
double *a;
ierr = VecGetArray(u, &a); CHKERRQ(ierr);
for (int i=0; i<m; i++) a[i] = exp(cos(i * PI / (m-1)));
for (int i=0; i<m1; i++) {
a[i] = exp(cos(i * PI / (m1-1)));
}
ierr = VecRestoreArray(u, &a); CHKERRQ(ierr);

// 2-D solution function
ierr = VecGetArray(u2, &a); CHKERRQ(ierr);
for (int i=0; i < m; i++) {
double x = (m==1) ? 0 : cos (i * PI / (m-1));
for (int j=0; j < n; j++) {
double y = (n==1) ? 0 : cos (j * PI / (n-1));
for (int k=0; k < p; k++) {
double z = (p==1) ? 0 : cos (k * PI / (p-1));
a[(i*n + j) * p + k] = exp(x) + exp(y) + exp(z);
}
}
}
ierr = VecRestoreArray(u2, &a); CHKERRQ(ierr);

ierr = MatMult(A, u, b); CHKERRQ(ierr);
ierr = MatMult(A2, u2, b2); CHKERRQ(ierr);
// ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
ierr = VecCopy(b, x); CHKERRQ(ierr);

Expand All @@ -60,10 +95,10 @@ int main(int argc,char **args)
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A\n",norm);CHKERRQ(ierr);

/* /\* ierr = VecView(x, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr); *\/ */
/* printf("\n"); */
/* ierr = VecView(b, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr); */
/* printf("\n"); */
/* ierr = VecView(u, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr); */
printf("\n");
ierr = VecView(b, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);
printf("\n");
ierr = VecView(b2, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);

/*
Free work space. All PETSc objects should be destroyed when they
Expand Down
149 changes: 149 additions & 0 deletions chebyshev.c
View file
@@ -1,4 +1,7 @@
#include "chebyshev.h"
#include <stdbool.h>

void perform_carry(ChebCtx *c, int *ind, int *offset, bool *done);

#undef __FUNCT__
#define __FUNCT__ "MatCreateChebD1"
Expand Down Expand Up @@ -80,3 +83,149 @@ PetscErrorCode ChebD1Destroy (Mat A) {
PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "MatCreateCheb"
PetscErrorCode MatCreateCheb(MPI_Comm comm, int rank, int tr, int *dim, unsigned flag,
Vec vx, Vec vy, Mat *A) {
PetscErrorCode ierr;
ChebCtx *c;
int n;
double *x, *y;

PetscFunctionBegin;
ierr = VecGetSize(vx, &n); CHKERRQ(ierr);
if (n < 2) SETERRQ1(PETSC_ERR_USER, "n = %d but must be >= 2", n);

ierr = PetscMalloc(sizeof(ChebCtx), &c); CHKERRQ(ierr);
ierr = PetscMalloc((rank - 1) * sizeof(fftw_iodim), &(c->dim)); CHKERRQ(ierr);
c->work = fftw_malloc(n * sizeof(double));
c->rank = rank;
c->tr = tr;

if (!(0 <= tr && tr < rank)) SETERRQ(PETSC_ERR_USER, "tdim out of range");
int stride = 1;
for (int r = rank-1, ri=rank-2; r >= 0; r--) {
fftw_iodim *iod;
if (r == tr) {
iod = &(c->tdim);
} else {
iod = &(c->dim[ri]);
ri--;
}
iod->n = dim[r];
iod->is = stride;
iod->os = stride;
stride *= dim[r];
}

if (n != stride) SETERRQ2(PETSC_ERR_USER, "dimensions do not agree: n = %d but stride = %d", n, stride);

ierr = VecGetArray(vx, &x); CHKERRQ(ierr);
ierr = VecGetArray(vy, &y); CHKERRQ(ierr);
const fftw_r2r_kind redft00 = FFTW_REDFT00, rodft00 = FFTW_RODFT00;
c->p_forward = fftw_plan_guru_r2r(1, &(c->tdim), c->rank-1, c->dim, x, c->work, &redft00, flag);
fftw_iodim tdim1 = c->tdim; tdim1.n -= 2; // we need a modified tdim to come back.
c->p_backward = fftw_plan_guru_r2r(1, &tdim1, c->rank-1, c->dim, c->work + tdim1.os, y + tdim1.is, &rodft00, flag);
ierr = VecRestoreArray(vx, &x); CHKERRQ(ierr);
ierr = VecRestoreArray(vy, &y); CHKERRQ(ierr);

ierr = MatCreateShell(comm, n, n, n, n, c, A); CHKERRQ(ierr);
ierr = MatShellSetOperation(*A, MATOP_MULT, (void(*)(void))ChebMult); CHKERRQ(ierr);
ierr = MatShellSetOperation(*A, MATOP_DESTROY, (void(*)(void))ChebDestroy); CHKERRQ(ierr);

PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "ChebMult"
PetscErrorCode ChebMult(Mat A, Vec vx, Vec vy) {
PetscErrorCode ierr;
double *x, *y;
ChebCtx *c;

PetscFunctionBegin;
ierr = MatShellGetContext(A, (void *)&c); CHKERRQ(ierr);
ierr = VecGetArray(vx, &x); CHKERRQ(ierr);
ierr = VecGetArray(vy, &y); CHKERRQ(ierr);

int n = c->tdim.n - 1; // Gauss-Lobatto-Chebyshev points are numbered from [0..n]
int ind[c->rank-1]; //ierr = PetscMalloc(c->rank * sizeof(int), &ind); CHKERRQ(ierr);

fftw_execute_r2r(c->p_forward, x, c->work);

ierr = PetscMemzero(ind, (c->rank - 1) * sizeof(int)); CHKERRQ(ierr);
int offset = 0;
for (bool done = false; !done; ) {
for (int i = 1; i < n; i++) {
int ix = offset + i * c->tdim.is;
c->work[ix] *= (double)i;
}
perform_carry(c, ind, &offset, &done);
}

fftw_execute_r2r(c->p_backward, c->work + c->tdim.os, y + c->tdim.is);

double N = (double)n;
double pin = PI / N;
ierr = PetscMemzero(ind, (c->rank - 1) * sizeof(int)); CHKERRQ(ierr);
offset = 0;
for (bool done = false; !done; ) {
int ix0 = offset;
int ixn = offset + n * c->tdim.is;
y[ix0] = 0.0;
y[ixn] = 0.0;
double s = 1.0;
for (int i = 1; i < n; i++) {
int ix = offset + i * c->tdim.is;
double I = (double)i;
y[ix] /= 2 * n * sqrt(1.0 - PetscSqr(cos(I * pin)));
y[ix0] += I * c->work[ix];
y[ixn] += s * I * c->work[ix];
s = -s;
}
y[ix0] = 0.5 * c->work[ixn] * N + y[ix0] / n;
y[ixn] = y[ixn] / N + 0.5 * s * N * c->work[ixn];
perform_carry(c, ind, &offset, &done);
}

ierr = VecRestoreArray(vx, &x); CHKERRQ(ierr);
ierr = VecRestoreArray(vy, &y); CHKERRQ(ierr);

PetscFunctionReturn(0);
}


void perform_carry(ChebCtx *c, int *ind, int *offset, bool *done) {
int carry;

*offset = 0;
carry = 1;
for (int i = c->rank-2; i >= 0; i--) {
ind[i] += carry;
if (ind[i] < c->dim[i].n) {
carry = 0;
} else {
ind[i] = 0;
carry = 1;
}
*offset += ind[i] * c->dim[i].is;
}
*done = (carry == 1);
}


#undef __FUNCT__
#define __FUNCT__ "ChebDestroy"
PetscErrorCode ChebDestroy (Mat A) {
PetscErrorCode ierr;
ChebCtx *c;

PetscFunctionBegin;
ierr = MatShellGetContext(A, (void *)&c); CHKERRQ(ierr);
fftw_destroy_plan(c->p_forward);
fftw_destroy_plan(c->p_backward);
fftw_free(c->work);
ierr = PetscFree(c->dim); CHKERRQ(ierr);
PetscFunctionReturn(0);
}

19 changes: 16 additions & 3 deletions chebyshev.h
View file
Expand Up @@ -2,9 +2,8 @@
#define CHEBYSHEV_H

// #include <math.h>
#include "fftw3.h"
#include "petscmat.h"
// #include "petscerror.h"
#include <fftw3.h>
#include <petscmat.h>

#define PI 3.14159265358979323846

Expand All @@ -14,8 +13,22 @@ typedef struct {
double *work;
} ChebD1Ctx;

typedef struct {
int rank, tr;
fftw_iodim tdim;
fftw_iodim *dim;
fftw_plan p_forward, p_backward;
double *work;
} ChebCtx;


PetscErrorCode MatCreateChebD1(MPI_Comm comm, Vec vx, Vec vy, unsigned flag, Mat *A);
PetscErrorCode ChebD1Mult(Mat A, Vec vx, Vec vy);
PetscErrorCode ChebD1Destroy (Mat A);

PetscErrorCode MatCreateCheb(MPI_Comm comm, int rank, int tr, int *dims, unsigned flag,
Vec vx, Vec vy, Mat *A);
PetscErrorCode ChebMult(Mat A, Vec vx, Vec vy);
PetscErrorCode ChebDestroy (Mat A);

#endif // CHEBYSHEV_H

0 comments on commit b2d95fb

Please sign in to comment.