cg.c
/* Least-squares solution using CG. The residual must be precomputed. */
Least-squares solution using CG. The residual must be precomputed.
Least - squares solution using CG . The residual must be precomputed .
real cgls (cplx *rn, cplx *sol, solveparm *slv,
		measdesc *src, measdesc *obs, real regparm) {
	long i, nelt = (long)fmaconf.numbases * (long)fmaconf.bspboxvol;
	int j, k;
	cplx *ifld, *pn, *scat, *adjcrt;
	real rnorm, pnorm, lrnorm, lpnorm, alpha, beta, rninc;
	ifld = malloc (3L * nelt * sizeof(cplx));
	pn = ifld + nelt;
	adjcrt = pn + nelt;
	scat = malloc (obs->count * sizeof(cplx));
	rnorm = 0;
	for (i = 0; i < nelt; ++i) {
		pn[i] = rn[i];
		sol[i] = 0.0;
		alpha = cabs (rn[i]);
		rnorm += alpha * alpha;
	}
	MPI_Allreduce (MPI_IN_PLACE, &rnorm, 1, MPIREAL, MPI_SUM, MPI_COMM_WORLD);
	pnorm = (rninc = rnorm);
	for (j = 0; j < slv->maxit; ++j) {
		memset (adjcrt, 0, nelt * sizeof(cplx));
		alpha = 0;
		for (i = 0; i < src->count; ++i) {
			buildrhs (ifld, src->locations + 3 * i, src->plane, NULL);
			bicgstab (ifld, ifld, 0, slv->maxit, slv->epscg, 1);
			frechet (pn, ifld, scat, obs, slv);
			frechadj (scat, ifld, adjcrt, obs, slv);
			for (k = 0; k < obs->count; ++k) {
				beta = cabs (scat[k]);
				alpha += beta * beta;
			}
		}
		alpha  = rnorm / (alpha + regparm * pnorm);
		beta = rnorm;
		for (i = 0, lrnorm = 0; i < nelt; ++i) {
			rn[i] -= alpha * (adjcrt[i] + regparm * pn[i]);
			rnorm = cabs (rn[i]);
			lrnorm += rnorm * rnorm;
		}
		MPI_Allreduce (&lrnorm, &rnorm, 1, MPIREAL, MPI_SUM, MPI_COMM_WORLD);
		beta = rnorm / beta;
		pnorm = lpnorm = 0;
		for (i = 0, lpnorm = 0; i < nelt; ++i) {
			sol[i] += alpha * pn[i];
			pn[i] = rn[i] + beta * pn[i];
			pnorm = cabs(pn[i]);
			lpnorm += pnorm * pnorm;
		}
		MPI_Allreduce (&lpnorm, &pnorm, 1, MPIREAL, MPI_SUM, MPI_COMM_WORLD);
		if (sqrt(rnorm / rninc) < slv->epscg) break;
	}
	free (ifld);
	free (scat);
	return sqrt(rnorm / rninc);
}
real cgls ( cplx * rn , cplx * sol , solveparm * slv , measdesc * src , measdesc * obs , real regparm ) { long i , nelt = ( long ) fmaconf . numbases * ( long ) fmaconf . bspboxvol ; int j , k ; cplx * ifld , * pn , * scat , * adjcrt ; real rnorm , pnorm , lrnorm , lpnorm , alpha , beta , rninc ; ifld = malloc ( 3L * nelt * sizeof ( cplx ) ) ; pn = ifld + nelt ; adjcrt = pn + nelt ; scat = malloc ( obs -> count * sizeof ( cplx ) ) ; rnorm = 0 ; for ( i = 0 ; i < nelt ; ++ i ) { pn [ i ] = rn [ i ] ; sol [ i ] = 0.0 ; alpha = cabs ( rn [ i ] ) ; rnorm += alpha * alpha ; } MPI_Allreduce ( MPI_IN_PLACE , & rnorm , 1 , MPIREAL , MPI_SUM , MPI_COMM_WORLD ) ; pnorm = ( rninc = rnorm ) ; for ( j = 0 ; j < slv -> maxit ; ++ j ) { memset ( adjcrt , 0 , nelt * sizeof ( cplx ) ) ; alpha = 0 ; for ( i = 0 ; i < src -> count ; ++ i ) { buildrhs ( ifld , src -> locations + 3 * i , src -> plane , NULL ) ; bicgstab ( ifld , ifld , 0 , slv -> maxit , slv -> epscg , 1 ) ; frechet ( pn , ifld , scat , obs , slv ) ; frechadj ( scat , ifld , adjcrt , obs , slv ) ; for ( k = 0 ; k < obs -> count ; ++ k ) { beta = cabs ( scat [ k ] ) ; alpha += beta * beta ; } } alpha = rnorm / ( alpha + regparm * pnorm ) ; beta = rnorm ; for ( i = 0 , lrnorm = 0 ; i < nelt ; ++ i ) { rn [ i ] -= alpha * ( adjcrt [ i ] + regparm * pn [ i ] ) ; rnorm = cabs ( rn [ i ] ) ; lrnorm += rnorm * rnorm ; } MPI_Allreduce ( & lrnorm , & rnorm , 1 , MPIREAL , MPI_SUM , MPI_COMM_WORLD ) ; beta = rnorm / beta ; pnorm = lpnorm = 0 ; for ( i = 0 , lpnorm = 0 ; i < nelt ; ++ i ) { sol [ i ] += alpha * pn [ i ] ; pn [ i ] = rn [ i ] + beta * pn [ i ] ; pnorm = cabs ( pn [ i ] ) ; lpnorm += pnorm * pnorm ; } MPI_Allreduce ( & lpnorm , & pnorm , 1 , MPIREAL , MPI_SUM , MPI_COMM_WORLD ) ; if ( sqrt ( rnorm / rninc ) < slv -> epscg ) break ; } free ( ifld ) ; free ( scat ) ; return sqrt ( rnorm / rninc ) ; }
Least-squares solution using CG.
Least - squares solution using CG .
/* Initialize some values. */ /* Reduce the local norms into a global norm. */ /* Build the incident field. */ /* Solve for the internal field. */ /* Compute the Frechet derivative. */ /* Compute the adjoint Frechet derivative. */ /* Update the local norm. */ /* Update the residual. */ /* Reduce the new residual norm. */ /* Update the solution and the search vector. */ /* Reduce the new search norm. */