#include <math.h>
#include <stdio.h>

/* #define TOEPLITZ */
#define ELLIPSE_5
/* #define ELLIPSE_7 */

#if defined(ELLIPSE_5)

#define SIZE_X           400
#define SIZE_Y           400
#define R               1.0 
#define PROB_SIZE       (SIZE_X * SIZE_Y)
#define MAX_NZERO       5

/*
  楕円型偏微分方程式の境界値問題

  -u_xx-u_yy+R*u_x = g(x,y)
  x:[0,1]  y:[0,1]

  厳密解 u(x,y)=1+xy) となるように
  g(x,y,z) を決める。

  ５点差分公式で解く
*/

void GetInfoProblem(int *size, int *nsize, int *max_nsize)
{
    *size=PROB_SIZE;
    *nsize=PROB_SIZE
          + 2 * (PROB_SIZE - SIZE_X)
          + 2 * (PROB_SIZE - SIZE_Y);
    *max_nsize=MAX_NZERO;
   }

void Problem(int i,int ja[],double a[],double *b,int *tsize)
{
    double inv_h = (SIZE_X + 1);
    double d0 =  4.0 * inv_h * inv_h;
    double d1 = -1.0 * inv_h * inv_h;
    double d2 = -(1.0 * inv_h * inv_h + 1.0 * inv_h * R / 2.0);
    double d3 = -(1.0 * inv_h * inv_h - 1.0 * inv_h * R / 2.0);
    double x,y,z;
    int k=0;
    if( i<0 || i>=PROB_SIZE) {
        printf("Error in Problem1");
        exit(1);
       }
    k=0;
    if (i >= SIZE_X ) {
       ja[k] = i - SIZE_X;
       a[k] = d1;
       ++k;
      }
    if (i % SIZE_X > 0) {
       ja[k] = i - 1;
       a[k] = d2;
       ++k;
      }
    ja[k] = i;
    a[k] = d0;
    ++k;
    if ((i + 1) % SIZE_X > 0) {
       ja[k] = i + 1;
       a[k] = d3;
       ++k;
      }
    if (i < PROB_SIZE - SIZE_X) {
       ja[k] = i + SIZE_X;
       a[k] = d1;
       ++k;
      }
    *tsize=k;
    x= (double)( i % SIZE_X + 1) / ( SIZE_X + 1 ) ;
    y= (double)( i / SIZE_Y + 1) / ( SIZE_Y + 1 ) ;

    *b = R * y;
    if( i<SIZE_X ) *b += inv_h * inv_h;
    if( i>=PROB_SIZE-SIZE_X ) *b += inv_h * inv_h * (1+x);
    if( i % SIZE_X == 0 ) *b += ( inv_h * inv_h + inv_h * R / 2.0);
    if( (i+1)% SIZE_X == 0 ) *b += ( inv_h * inv_h - inv_h * R / 2.0)*(1+y);
}
#endif


#if defined(ELLIPSE_7)

#define SIZE_X          80
#define SIZE_Y          80
#define SIZE_Z          80
#define R               100.0
#define PROB_SIZE       (SIZE_X * SIZE_Y * SIZE_Z)
#define SIZE_XY         (SIZE_X * SIZE_Y)
#define MAX_NZERO       7
#define PI              3.141592653589793

/*
  楕円型偏微分方程式の境界値問題

  -u_xx-u_yy-u_zz+R*u_x = g(x,y,z)
  x:[0,1]  y:[0,1]  z:[0,1]

  厳密解 u(x,y,z)=exp(x*y*z)*sin(PI*x)*sin(PI*y)*sin(PI*z) となるように
  g(x,y,z) を決める。

  ７点差分公式で解く
*/

void GetInfoProblem(int *size, int *nsize, int *max_nsize)
{
    *size=PROB_SIZE;
    *nsize=PROB_SIZE
          + 2 * (PROB_SIZE - SIZE_XY)
	  + 2 * (PROB_SIZE - SIZE_Y * SIZE_Z)
	  + 2 * (PROB_SIZE - SIZE_Z * SIZE_X);
    *max_nsize=MAX_NZERO;
}

void Problem(int i,int ja[],double a[],double *b,int *tsize)
{
    double inv_h = (SIZE_X + 1); 
    double d0 =  6.0 * inv_h * inv_h;
    double d1 = -1.0 * inv_h * inv_h;
    double d2 = -(1.0 * inv_h * inv_h + 1.0 * inv_h * R / 2.0);
    double d3 = -(1.0 * inv_h * inv_h - 1.0 * inv_h * R / 2.0);
    double x,y,z;
    int k=0;

    if( i<0 || i>=PROB_SIZE) {
        printf("Error in Problem1");
        exit(1);
    }
    k=0;
    if (i >= SIZE_XY) {
       ja[k] = i - SIZE_XY;
       a[k] = d1;
       ++k;
    }
    if(i % SIZE_XY >= SIZE_X){
       ja[k] = i - SIZE_X;
       a[k] = d1;
       ++k;
    }
    if (i % SIZE_XY % SIZE_X > 0) {
       ja[k] = i - 1;
       a[k] = d2;
       ++k;
    }
    ja[k] = i;
    a[k] = d0;
    ++k;
    
    if ((i + 1) % SIZE_XY % SIZE_X > 0) {
       ja[k] = i + 1;
       a[k] = d3;
       ++k;
    }
    if (i % SIZE_XY < SIZE_XY - SIZE_X) {
       ja[k] = i + SIZE_X;
       a[k] = d1;
       ++k;
    }
    if (i < PROB_SIZE - SIZE_XY) {
       ja[k] = i + SIZE_XY;
       a[k] = d1;
       ++k;
    }
    *tsize=k;
    x= (double)( i % SIZE_XY % SIZE_X + 1) / ( SIZE_X + 1 ) ;
    y= (double)( i % SIZE_XY / SIZE_Y + 1) / ( SIZE_Y + 1 ) ;
    z= (double)( i / SIZE_XY + 1) / ( SIZE_Z + 1 ) ;

    *b = -y*y*z*z*exp(x*y*z)*sin(PI*x)*sin(PI*y)*sin(PI*z)
       - 2*PI*y*z*exp(x*y*z)*cos(PI*x)*sin(PI*y)*sin(PI*z)
       - x*x*z*z*exp(x*y*z)*sin(PI*x)*sin(PI*y)*sin(PI*z)
       - 2*PI*x*z*exp(x*y*z)*sin(PI*x)*cos(PI*y)*sin(PI*z)
       - x*x*y*y*exp(x*y*z)*sin(PI*x)*sin(PI*y)*sin(PI*z)
       - 2*PI*x*y*exp(x*y*z)*sin(PI*x)*sin(PI*y)*cos(PI*z)
       + 3*PI*PI*exp(x*y*z)*sin(PI*x)*sin(PI*y)*sin(PI*z)
       + R*(y*z*exp(x*y*z)*sin(PI*x)*sin(PI*y)*sin(PI*z)
            +PI*exp(x*y*z)*cos(PI*x)*sin(PI*y)*sin(PI*z));
}
#endif

#if defined(TOEPLITZ)

#define SIZE           4000000
#define R              1.0
#define PROB_SIZE      SIZE
#define MAX_NZERO      3

/*
  Toeplitz行列
  2 1 0
  0 2 1
  r 0 2 1
    r 0 2 1
     ・・・
        r 0 2
*/

void GetInfoProblem(int *size, int *nsize, int *max_nsize)
{
    *size=PROB_SIZE;
    *nsize=PROB_SIZE*3-3;
    *max_nsize=MAX_NZERO;
}

void Problem(int i,int ja[],double a[],double *b,int *tsize)
{
    int k=0;

    if( i<0 || i>=PROB_SIZE) {
        printf("Error in Problem1");
        exit(1);
    }
    k=0;
    if (i >= 2) {
       ja[k] = i - 2;
       a[k] = R;
       ++k;
    }
    ja[k] = i;
    a[k] = 2.0;
    ++k;
    if (i < SIZE-1) {
       ja[k] = i + 1;
       a[k] = 1.0;
       ++k;
    }
    *tsize=k;
    *b = 1.0;
}
#endif