frechet.c
/* Estimate the square of the spectral radius of the Frechet derivative
 * using power iteration on the product of the derivative and its adjoint. */
Estimate the square of the spectral radius of the Frechet derivative
using power iteration on the product of the derivative and its adjoint.
Estimate the square of the spectral radius of the Frechet derivative using power iteration on the product of the derivative and its adjoint .
real specrad (int maxit, solveparm *slv, measdesc *src, measdesc *obs) {
	cplx *ifld, *adjcrt, *scat, *pn, mu = 1.0, newmu;
	long k, nelt = (long)fmaconf.numbases * (long)fmaconf.bspboxvol,
	     randmax = (2L << 31) - 1;
	real nrm, tmp;
	int i, j;
	ifld = malloc (3L * nelt * sizeof(cplx));
	adjcrt = ifld + nelt;
	pn = adjcrt + nelt;
	scat = malloc (obs->count * sizeof(cplx));
	for (k = 0, nrm = 0.0; k < nelt; ++k) {
		pn[k] = ((real)random() + I * (real)random()) / (real)randmax;
		tmp = cabs(pn[k]);
		nrm += tmp * tmp;
	}
	MPI_Allreduce (MPI_IN_PLACE, &nrm, 1, MPIREAL, MPI_SUM, MPI_COMM_WORLD);
	nrm = sqrt(nrm);
	for (k = 0; k < nelt; ++k) pn[k] /= nrm;
	for (j = 0; j < maxit; ++j) {
		memset (adjcrt, 0, nelt * sizeof(cplx));
		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, nrm = newmu = 0.0; k < nelt; ++k) {
			newmu += conj(pn[k]) * adjcrt[k];
			tmp = cabs(adjcrt[k]);
			nrm += tmp * tmp;
		}
		MPI_Allreduce (MPI_IN_PLACE, &nrm, 1, MPIREAL, MPI_SUM, MPI_COMM_WORLD);
		MPI_Allreduce (MPI_IN_PLACE, &newmu, 2, MPIREAL, MPI_SUM, MPI_COMM_WORLD);
		nrm = sqrt(nrm);
		if (cabs((newmu - mu) / mu) < slv->epscg) break;
		for (k = 0; k < nelt; ++k)
			pn[k] = adjcrt[k] / nrm;
		mu = newmu;
	}
	free (ifld);
	free (scat);
	return creal(mu);
}
real specrad ( int maxit , solveparm * slv , measdesc * src , measdesc * obs ) { cplx * ifld , * adjcrt , * scat , * pn , mu = 1.0 , newmu ; long k , nelt = ( long ) fmaconf . numbases * ( long ) fmaconf . bspboxvol , randmax = ( 2L << 31 ) - 1 ; real nrm , tmp ; int i , j ; ifld = malloc ( 3L * nelt * sizeof ( cplx ) ) ; adjcrt = ifld + nelt ; pn = adjcrt + nelt ; scat = malloc ( obs -> count * sizeof ( cplx ) ) ; for ( k = 0 , nrm = 0.0 ; k < nelt ; ++ k ) { pn [ k ] = ( ( real ) random ( ) + I * ( real ) random ( ) ) / ( real ) randmax ; tmp = cabs ( pn [ k ] ) ; nrm += tmp * tmp ; } MPI_Allreduce ( MPI_IN_PLACE , & nrm , 1 , MPIREAL , MPI_SUM , MPI_COMM_WORLD ) ; nrm = sqrt ( nrm ) ; for ( k = 0 ; k < nelt ; ++ k ) pn [ k ] /= nrm ; for ( j = 0 ; j < maxit ; ++ j ) { memset ( adjcrt , 0 , nelt * sizeof ( cplx ) ) ; 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 , nrm = newmu = 0.0 ; k < nelt ; ++ k ) { newmu += conj ( pn [ k ] ) * adjcrt [ k ] ; tmp = cabs ( adjcrt [ k ] ) ; nrm += tmp * tmp ; } MPI_Allreduce ( MPI_IN_PLACE , & nrm , 1 , MPIREAL , MPI_SUM , MPI_COMM_WORLD ) ; MPI_Allreduce ( MPI_IN_PLACE , & newmu , 2 , MPIREAL , MPI_SUM , MPI_COMM_WORLD ) ; nrm = sqrt ( nrm ) ; if ( cabs ( ( newmu - mu ) / mu ) < slv -> epscg ) break ; for ( k = 0 ; k < nelt ; ++ k ) pn [ k ] = adjcrt [ k ] / nrm ; mu = newmu ; } free ( ifld ) ; free ( scat ) ; return creal ( mu ) ; }
Estimate the square of the spectral radius of the Frechet derivative
 using power iteration on the product of the derivative and its adjoint.
Estimate the square of the spectral radius of the Frechet derivative using power iteration on the product of the derivative and its adjoint .
/* Build a random initial vector. */ /* Normalize the initial vector. */ /* Loop through power iterations. */ /* Compute the product (F*)(F)p. */ /* Build the incident field. */ /* Solve for the internal field. */ /* Compute the Frechet derivative. */ /* Compute the adjoint Frechet derivative. */ /* Update the norm of the solution and find the approximate
		 * singular value. Note that the norm of the test vector is 1. */ /* Break if the updated singular value hasn't changed much. */ /* Normalize the solution for the next iteration. */ /* Copy the new singular value. */