Implementation of alternating direction implicit method

I have written a program that implements the ADI method and Crank-Nicolson method for solving Schrodinger equations.

The program is working, but it takes a very long time to run. I am looking for tips on how to improve the performance of the program.

```
#define _USE_MATH_DEFINES
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include

#define PI M_PI
#define h_bar 1.0
#define M 1.0
#define L .5
#define sig .05
#define grid_point 100
#define K1 1.0E6
#define K2 1.0E6
#define omega1 std::sqrt(K1/M)
#define omega2 std::sqrt(K2/M)

typedef std::complex compx;
#define I compx(0.,1.)
#define one compx(1.,0.)
#define two compx(2.,0.)

class wave_function {
public:
wave_function();
wave_function(bool);
double dt = 4.2E-06;
double dx = L / grid_point;
std::vector> value;
void solve_triag();
double potential(int, int);
double real_space(int);
double density_x1[grid_point];
double density_x2[grid_point];
compx sec_space_deriv(char, int, int);
void normalize();
void rho();
void alpha_beta_solver(wave_function &, compx, std::vector &mid, compx, std::vector &R, int, char);
compx a = -h_bar / (twoM), b = -h_bar / (twoM);
double r = dt / (dx*dx);
compx A = (Ira / two), B = (Irb / two), C = I*dt / two;
};

wave_function::wave_function() {
value.resize(grid_point);
for (int l = 0; l mid1, mid2, R1, R2;
wave_function tmp;

mid1.resize(grid_point);
R1.resize(grid_point);
for (int x2 = 0; x2 &mid, compx post_mid, std::vector &R, int coor, char x_or_y) {

compx x_N, R_N;// new_mid[grid_point], new_R[grid_point];
std::vector alpha(grid_point - 2);
std::vector beta1(grid_point - 1);//for x_1
std::vector beta2(grid_point - 1);//for x_2
std::vector x_1(grid_point), x_2(grid_point);

alpha[0] = post_mid / mid[0];
beta1[0] = R[0] / mid[0];
beta2[0] = -b4_mid / mid[0];

//Forward run
for (int l = 1; l = 0; l--) {

First, you should fix these compiler warnings:

157094.cpp: In constructor ‘wave_function::wave_function()’:
157094.cpp:53:1: warning: ‘wave_function::value’ should be initialized in the member initialization list [-Weffc++]
wave_function::wave_function() {
^~~~~~~~~~~~~
157094.cpp: In constructor ‘wave_function::wave_function(bool)’:
157094.cpp:63:1: warning: ‘wave_function::value’ should be initialized in the member initialization list [-Weffc++]
wave_function::wave_function(bool init) {
^~~~~~~~~~~~~
157094.cpp: In member function ‘void wave_function::alpha_beta_solver(wave_function&, compx, std::vector >&, compx, std::vector >&, int, char)’:
157094.cpp:119:54: warning: unused parameter ‘New’ [-Wunused-parameter]
void wave_function::alpha_beta_solver(wave_function &New, compx b4_mid, std::vector &mid, compx post_mid, std::vector &R, int coor, char x_or_y) {
^~~
157094.cpp: In function ‘int main()’:
157094.cpp:249:10: warning: unused variable ‘check’ [-Wunused-variable]
double check = 0.0;
^~~~~
157094.cpp: In member function ‘compx wave_function::sec_space_deriv(char, int, int)’:
157094.cpp:205:1: warning: control reaches end of non-void function [-Wreturn-type]
}
^


I'd advise that you use real, strongly-typed constants in preference to preprocessor substitution:

const double hbar{1};
const double L{0.5};
// etc.

const compx I{0,1};
// No need for complex values of 1 or 2, as doubles will promote okay

The following code doesn't do what you think it does:

wave_function::wave_function(bool init) {
  if (init) {
    // (snip)
  }
  else if (!init) {
    wave_function();
  }
}

In the else branch (with a redundant condition, but your compiler can probably eliminate that), you construct a new wave_function object and then discard it. wave_function(); here will not invoke the default constructor.

I think you want to turn this the other way around:

wave_function::wave_function()
    : wave_function(false)
{}

wave_function::wave_function(bool init) {
  if (init) {
      // "true" setup
  }
  else {
      // "false" setup
  }
}

The population of the initial values is hard to read, and repeats a lot of calculations. These can be avoided with the judicious use of temporary variables for the intermediate results. I believe the following is equivalent to your initialization loop:

for (int l = 0; l < grid_point; l++) {
        value[l].resize(grid_point);
        for (int m = 0; m < grid_point; m++) {
            auto const rl = real_space(l);
            auto const rm = real_space(m);

            auto const k1 = M * omega1 / PI / h_bar;
            auto const k2 = M * omega2 / PI / h_bar;

            auto const e1 = std::exp(-M*omega1*rl*rl/h_bar/2);
            auto const e2 = std::exp(-M*omega2*rm*rm/h_bar/2);

            auto const psi00 = std::pow(k1*k2, .25) * e1 * e2;
            auto const psi10 = psi00 * std::sqrt(2*M*omega1/h_bar) * rl;
            value[l][m] = psi00+psi10;
        }
    }

I've used the identity pow(x,n)pow(y,n) == pow(xy,n) in a couple of places above, to reduce the number of calls.

Throughout your code, there are uses of std::pow with constant exponents. In general, you can expect x*x to be faster than std::pow(x, 2) and std::sqrt(x) to be faster than std::pow(x, 0.5) - and both are likely to be more accurate into the bargain.

(This is an incomplete answer - Real Life called me away. I was intending to return with additional review, but Edward has written an answer much better than this one was ever going to be.)

StackExchange Code Review Q#157094, answer score: 2