mlfma.c
/* Construct an ACA approximation to the far-field signature matrix. */
Construct an ACA approximation to the far-field signature matrix.
Construct an ACA approximation to the far - field signature matrix .
static int acabuild (cplx **mats, real k0, real tol,
		int ntheta, int nphi, real dx, int bpd) {
	int maxrank, *irow, *icol, i, j, k,
	    nsamp = (ntheta - 2) * nphi + 2, nelt = bpd * bpd * bpd;
	cplx *u, *v, *uptr, *vptr, *row, *col, dpr, dpc;
	real dist[3], s[3], err = 0, tolsq;
	tolsq = tol * tol * 1e-2;
	maxrank = MIN(nelt, nsamp);
	irow = malloc((2 * maxrank + 1) * sizeof(int));
	icol = irow + maxrank + 1;
	u = malloc(maxrank * (nelt + nsamp) * sizeof(cplx));
	v = u + nsamp * maxrank;
	irow[0] = 0;
	for (i = 0, row = v, col = u; i < maxrank; ++i, row += nelt, col += nsamp) {
		sampcoords (s, irow[i], ntheta, nphi);
		farmatrow (row, k0, s, dx, bpd);
#pragma omp parallel for default(shared) private(vptr, uptr, j, k)
		for (j = 0; j < nelt; ++j) {
			for (k = 0; k < i; ++k, uptr += nsamp, vptr += nelt) {
				uptr = u + k * nsamp;
				vptr = v + k * nelt;
				row[j] -= vptr[j] * uptr[irow[i]];
			}
		}
		icol[i] = maxind (row, nelt, icol, i);
		dpr = row[icol[i]];
#pragma omp parallel for default(shared) private(j)
		for (j = 0; j < nelt; ++j) row[j] /= dpr;
		cellcoords (dist, icol[i], bpd, dx);
		farmatcol (col, k0, dist, ntheta, nphi);
#pragma omp parallel for default(shared) private(vptr, uptr, j, k)
		for (j = 0; j < nsamp; ++j) {
			for (k = 0; k < i; ++k) {
				vptr = v + k * nelt;
				uptr = u + k * nsamp;
				col[j] -= vptr[icol[i]] * uptr[j];
			}
		}
		for (k = 0, uptr = u, vptr = v; k < i; ++k) {
			DOTU_SUB (nelt, vptr, 1, row, 1, &dpr);
			DOTU_SUB (nsamp, uptr, 1, col, 1, &dpc);
			uptr += nsamp;
			vptr += nelt;
			err += 2.0 * cabs(dpr) * cabs(dpc);
		}
		DOTC_SUB (nelt, row, 1, row, 1, &dpr);
		DOTC_SUB (nsamp, col, 1, col, 1, &dpc);
		err += creal(dpr) * creal(dpc);
		if (creal(dpr) * creal(dpc) <= tolsq * err) break;
		irow[i + 1] = maxind (col, nsamp, irow, i + 1);
	}
	maxrank = i;
	for (i = 0, k = maxrank * nelt; i < k; ++i) v[i] = conj(v[i]);
	maxrank = recaca (u, v, nsamp, nelt, maxrank, tol);
	*mats = malloc (maxrank * (nelt + nsamp) * sizeof(cplx));
	memcpy (*mats, u, maxrank * nsamp * sizeof(cplx));
	memcpy (*mats + maxrank * nsamp, v, maxrank * nelt * sizeof(cplx));
	free (irow);
	free (u);
	return maxrank;
}
static int acabuild ( cplx * * mats , real k0 , real tol , int ntheta , int nphi , real dx , int bpd ) { int maxrank , * irow , * icol , i , j , k , nsamp = ( ntheta - 2 ) * nphi + 2 , nelt = bpd * bpd * bpd ; cplx * u , * v , * uptr , * vptr , * row , * col , dpr , dpc ; real dist [ 3 ] , s [ 3 ] , err = 0 , tolsq ; tolsq = tol * tol * 1e-2 ; maxrank = MIN ( nelt , nsamp ) ; irow = malloc ( ( 2 * maxrank + 1 ) * sizeof ( int ) ) ; icol = irow + maxrank + 1 ; u = malloc ( maxrank * ( nelt + nsamp ) * sizeof ( cplx ) ) ; v = u + nsamp * maxrank ; irow [ 0 ] = 0 ; for ( i = 0 , row = v , col = u ; i < maxrank ; ++ i , row += nelt , col += nsamp ) { sampcoords ( s , irow [ i ] , ntheta , nphi ) ; farmatrow ( row , k0 , s , dx , bpd ) ; #pragma  omp parallel for default(shared) private(vptr, uptr, j, k) for ( j = 0 ; j < nelt ; ++ j ) { for ( k = 0 ; k < i ; ++ k , uptr += nsamp , vptr += nelt ) { uptr = u + k * nsamp ; vptr = v + k * nelt ; row [ j ] -= vptr [ j ] * uptr [ irow [ i ] ] ; } } icol [ i ] = maxind ( row , nelt , icol , i ) ; dpr = row [ icol [ i ] ] ; #pragma  omp parallel for default(shared) private(j) for ( j = 0 ; j < nelt ; ++ j ) row [ j ] /= dpr ; cellcoords ( dist , icol [ i ] , bpd , dx ) ; farmatcol ( col , k0 , dist , ntheta , nphi ) ; #pragma  omp parallel for default(shared) private(vptr, uptr, j, k) for ( j = 0 ; j < nsamp ; ++ j ) { for ( k = 0 ; k < i ; ++ k ) { vptr = v + k * nelt ; uptr = u + k * nsamp ; col [ j ] -= vptr [ icol [ i ] ] * uptr [ j ] ; } } for ( k = 0 , uptr = u , vptr = v ; k < i ; ++ k ) { DOTU_SUB ( nelt , vptr , 1 , row , 1 , & dpr ) ; DOTU_SUB ( nsamp , uptr , 1 , col , 1 , & dpc ) ; uptr += nsamp ; vptr += nelt ; err += 2.0 * cabs ( dpr ) * cabs ( dpc ) ; } DOTC_SUB ( nelt , row , 1 , row , 1 , & dpr ) ; DOTC_SUB ( nsamp , col , 1 , col , 1 , & dpc ) ; err += creal ( dpr ) * creal ( dpc ) ; if ( creal ( dpr ) * creal ( dpc ) <= tolsq * err ) break ; irow [ i + 1 ] = maxind ( col , nsamp , irow , i + 1 ) ; } maxrank = i ; for ( i = 0 , k = maxrank * nelt ; i < k ; ++ i ) v [ i ] = conj ( v [ i ] ) ; maxrank = recaca ( u , v , nsamp , nelt , maxrank , tol ) ; * mats = malloc ( maxrank * ( nelt + nsamp ) * sizeof ( cplx ) ) ; memcpy ( * mats , u , maxrank * nsamp * sizeof ( cplx ) ) ; memcpy ( * mats + maxrank * nsamp , v , maxrank * nelt * sizeof ( cplx ) ) ; free ( irow ) ; free ( u ) ; return maxrank ; }
Construct an ACA approximation to the far-field signature matrix.
Construct an ACA approximation to the far - field signature matrix .
/* Compute the square of the tolerance and decrease by two orders
	 * of magnitude for recompression. */ /* Allocate the row and column index arrays. The extra element should
	 * only be stored (but never recalled) in the limiting case when the
	 * rank cannot be reduced. */ /* Allocate the workspace for the row and column matrices. */ /* Start with the first row of the matrix. */ /* Find the coordinate of the observer element. */ /* Fill the row for the selected observer element. */ /* Subtract existing contributions from earlier ranks. */ /* Find the index and value of the maximum value. */ /* Scale the row. */ /* Find the coordinates of the source cell. */ /* Construct the radiation pattern of the source cell. */ /* Subtract existing contributions from earlier ranks. */ /* Update the error approximation. */ /* Find the next row index. */ /* Conjugate the matrix v. */ /* Recompress (in place) the matrices with the actual tolerance. */ /* Allocate the final matrix storage. */ /* Copy the colum matrix in first, then the row matrix. */ /* Free the work arrays. */