| 1) https://github.com/bcrestel/parallelmultigrid 2) void v_cycle(double*& array3d, double h, int N, double* function_3d, int v1, int v2) { std::cout << "v_cycle"<< std::endl; // base case if(N == 1){ smoother7p_3dfun(array3d, N, function_3d );} else{ // (1) pre-smooth (v1) times: for(int i = 0; i < v1; i++) smoother7p_3dfun(array3d, N, function_3d ); // (2) compute residual vector: double* mat_vec = matvecAv7p(array3d, N); double *residual = vecSub_fun3d(mat_vec, function_3d, N); delete[] mat_vec; // (3) restrict residual: restrict_3d_full_weighting(residual, N); // N will become floor(N/2) // (4) make the initial guess for subsequent calls = 0 (as we are solving for the error !) double* e_2h = new double[N*N*N]; #pragma omp parallel for for(int i = 0; i < (N*N*N); i++){ e_2h[i] = 0;} // (5) recurse: v_cycle(e_2h, 2*h, N, residual, v1, v2); delete[] residual; // (6) prolongate: prolongation3d(e_2h, N); // N will become N = 2*N+1 double* addition = vecAdd(array3d, e_2h, N); delete[] array3d; array3d = addition; // (7) post-smooth (v2) times .. for(int i = 0; i < v2; i++) smoother7p_3dfun(array3d, N, function_3d ); printf("v_cycle\n"); } } 3)如上所示，全陷入死循环，就没有搞了。 |
说点什么...