view release on metacpan or  search on metacpan

src/Source/FreeImageToolkit/MultigridPoissonSolver.cpp  view on Meta::CPAN

	const float *uc_bits = (float*)FreeImage_GetBits(UC);
	
	// do elements that are copies
	{
		const int nc = nf/2 + 1;

		float *uf_scan = uf_bits;
		const float *uc_scan = uc_bits;		
		for (row_uc = 0; row_uc < nc; row_uc++) {
			for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
				// calculate UF(2*row_uc, col_uf) = UC(row_uc, col_uc);
				uf_scan[col_uf] = uc_scan[col_uc];
			}
			uc_scan += uc_pitch;
			uf_scan += 2 * uf_pitch;
		}
	}
	// do odd-numbered columns, interpolating vertically
	{		
		for(row_uf = 1; row_uf < nf-1; row_uf += 2) {
			float *uf_scan = uf_bits + row_uf * uf_pitch;
			for (col_uf = 0; col_uf < nf; col_uf += 2) {
				// calculate UF(row_uf, col_uf) = 0.5 * ( UF(row_uf+1, col_uf) + UF(row_uf-1, col_uf) )
				uf_scan[col_uf] = 0.5F * ( *(uf_scan + uf_pitch + col_uf) + *(uf_scan - uf_pitch + col_uf) );
			}
		}
	}
	// do even-numbered columns, interpolating horizontally
	{
		float *uf_scan = uf_bits;
		for(row_uf = 0; row_uf < nf; row_uf++) {
			for (col_uf = 1; col_uf < nf-1; col_uf += 2) {
				// calculate UF(row_uf, col_uf) = 0.5 * ( UF(row_uf, col_uf+1) + UF(row_uf, col_uf-1) )
				uf_scan[col_uf] = 0.5F * ( uf_scan[col_uf + 1] + uf_scan[col_uf - 1] );
			}
			uf_scan += uf_pitch;
		}
	}
}

/**
Red-black Gauss-Seidel relaxation for model problem. Updates the current value of the solution
u[0..n-1][0..n-1], using the right-hand side function rhs[0..n-1][0..n-1].
*/
static void fmg_relaxation(FIBITMAP *U, FIBITMAP *RHS, int n) {
	int row, col, ipass, isw, jsw;
	const float h = 1.0F / (n - 1);
	const float h2 = h*h;

	const int u_pitch  = FreeImage_GetPitch(U) / sizeof(float);
	const int rhs_pitch  = FreeImage_GetPitch(RHS) / sizeof(float);
	
	float *u_bits = (float*)FreeImage_GetBits(U);
	const float *rhs_bits = (float*)FreeImage_GetBits(RHS);

	for (ipass = 0, jsw = 1; ipass < 2; ipass++, jsw = 3-jsw) { // Red and black sweeps
		float *u_scan = u_bits + u_pitch;
		const float *rhs_scan = rhs_bits + rhs_pitch;
		for (row = 1, isw = jsw; row < n-1; row++, isw = 3-isw) {
			for (col = isw; col < n-1; col += 2) { 
				// Gauss-Seidel formula
				// calculate U(row, col) = 
				// 0.25 * [ U(row+1, col) + U(row-1, col) + U(row, col+1) + U(row, col-1) - h2 * RHS(row, col) ]		 
				float *u_center = u_scan + col;
				const float *rhs_center = rhs_scan + col;
				*u_center = *(u_center + u_pitch) + *(u_center - u_pitch) + *(u_center + 1) + *(u_center - 1);
				*u_center -= h2 * *rhs_center;
				*u_center *= 0.25F;
			}
			u_scan += u_pitch;
			rhs_scan += rhs_pitch;
		}
	}
}

/**
Returns minus the residual for the model problem. Input quantities are u[0..n-1][0..n-1] and
rhs[0..n-1][0..n-1], while res[0..n-1][0..n-1] is returned.
*/
static void fmg_residual(FIBITMAP *RES, FIBITMAP *U, FIBITMAP *RHS, int n) {
	int row, col;

	const float h = 1.0F / (n-1);	
	const float h2i = 1.0F / (h*h);

	const int res_pitch  = FreeImage_GetPitch(RES) / sizeof(float);
	const int u_pitch  = FreeImage_GetPitch(U) / sizeof(float);
	const int rhs_pitch  = FreeImage_GetPitch(RHS) / sizeof(float);
	
	float *res_bits = (float*)FreeImage_GetBits(RES);
	const float *u_bits = (float*)FreeImage_GetBits(U);
	const float *rhs_bits = (float*)FreeImage_GetBits(RHS);

	// interior points
	{
		float *res_scan = res_bits + res_pitch;
		const float *u_scan = u_bits + u_pitch;
		const float *rhs_scan = rhs_bits + rhs_pitch;
		for (row = 1; row < n-1; row++) {
			for (col = 1; col < n-1; col++) {
				// calculate RES(row, col) = 
				// -h2i * [ U(row+1, col) + U(row-1, col) + U(row, col+1) + U(row, col-1) - 4 * U(row, col) ] + RHS(row, col);
				float *res_center = res_scan + col;
				const float *u_center = u_scan + col;
				const float *rhs_center = rhs_scan + col;
				*res_center = *(u_center + u_pitch) + *(u_center - u_pitch) + *(u_center + 1) + *(u_center - 1) - 4 * *u_center;
				*res_center *= -h2i;
				*res_center += *rhs_center;
			}
			res_scan += res_pitch;
			u_scan += u_pitch;
			rhs_scan += rhs_pitch;
		}
	}

	// boundary points
	{
		memset(FreeImage_GetScanLine(RES, 0), 0, FreeImage_GetPitch(RES));
		memset(FreeImage_GetScanLine(RES, n-1), 0, FreeImage_GetPitch(RES));
		float *left = res_bits;
		float *right = res_bits + (n-1);