matrix_algebra.cpp

矩量法仿真电磁辐射和散射的源代码（c++）

		}
	}
	
#ifdef NEC_MATRIX_CHECK
	cout << "atlas_solved = ";
	to_octave(a_in,n,ndim);

	cout << "atlas_ip = ";
	to_octave(ip,n);
#endif
} 
#endif /*  ATLAS */


/*-----------------------------------------------------------------------*/

/*	factrs

	For symmetric structure, transforms submatricies to form
	matricies of the symmetric modes and calls routine to LU decompose
	matricies.
	
	If no symmetry [nrow = np], the routine is called to LU decompose the
	complete matrix.
*/
void factrs(nec_output_file& s_output,  int np, int nrow, complex_array& a, int_array& ip )
{
	DEBUG_TRACE("factrs(" << np << "," << nrow << ")");
	if (nrow == np) // no symmetry
	{
		lu_decompose(s_output,  np, a, ip, nrow );
		return;
	}
	
	int num_symmetric_modes = nrow / np;
	DEBUG_TRACE("\tnum_symmetric_modes = " << num_symmetric_modes);
	
	for (int mode = 0; mode < num_symmetric_modes; mode++ )
	{
		int mode_offset = mode * np;
		
		complex_array a_temp = a.sub_array(mode_offset);
		int_array ip_temp = ip.sub_array(mode_offset);
		
		lu_decompose(s_output,  np, a_temp, ip_temp, nrow );
	}
}

/*-----------------------------------------------------------------------*/

/*! \brief
	Subroutine to solve the matrix equation lu*x=b where l is a unit
	lower triangular matrix and u is an upper triangular matrix both
	of which are stored in a.  the rhs vector b is input and the
	solution is returned through vector b.   (matrix transposed)
      
	  COMPLEX*16 A,B,Y,SUM
      INTEGER PI
      COMMON /SCRATM/ Y(2*MAXSEG)
      DIMENSION A(NDIM,NDIM), IP(NDIM), B(NDIM)
C
C     FORWARD SUBSTITUTION
C
      DO 3 I=1,N
      PI=IP(I)
      Y(I)=B(PI)
      B(PI)=B(I)
      IP1=I+1
      IF (IP1.GT.N) GO TO 2
      DO 1 J=IP1,N
      B(J)=B(J)-A(J,I)*Y(I)
1     CONTINUE
2     CONTINUE
3     CONTINUE
C
C     BACKWARD SUBSTITUTION
C
      DO 6 K=1,N
      I=N-K+1
      SUM=(0.,0.)
      IP1=I+1
      IF (IP1.GT.N) GO TO 5
      DO 4 J=IP1,N
      SUM=SUM+A(I,J)*B(J)
4     CONTINUE
5     CONTINUE
      B(I)=(Y(I)-SUM)/A(I,I)
6     CONTINUE
      RETURN
      END
	
*/
void solve( int n, complex_array& a, int_array& ip,
    complex_array& b, int ndim )
{
	DEBUG_TRACE("solve(" << n << "," << ndim << ")");
/*	
	We should use zgetrs from LAPACK to solve this.
*/	
	complex_array y(n);
	
	/* forward substitution */
	for (int i = 0; i < n; i++ )
	{
		int pivot_index = ip[i]-1;
		y[i]= b[pivot_index];
		b[pivot_index]= b[i];
		int ip1= i+1;
		
		int i_offset = i*ndim;
		for (int j = ip1; j < n; j++ )
			b[j] -= a[j+i_offset] * y[i];
	}
	
	/* backward substitution */
	for (int k = 0; k < n; k++ )
	{
		int i= n-k-1;
		nec_complex sum(cplx_00());
		int ip1= i+1;
		
		for (int j = ip1; j < n; j++ )
			sum += a[i+j*ndim]* b[j];
		
		b[i] = (y[i]- sum) / a[i+i*ndim];
	}
}



/*
	Subroutine solves, for symmetric structures, handles the
	transformation of the right hand side vector and solution
	of the matrix eq.
*/
void solves(complex_array& a, int_array& ip, complex_array& b, int neq,
	int nrh, int np, int n, int mp, int m, int nop, 
	complex_array& symmetry_array
)
{
	DEBUG_TRACE("solves(" << neq << "," << nrh << "," << np << "," << n << ")");
	int  ic, i, kk, ia, ib, j, k;

	nec_complex  sum;
	
	/* Allocate to scratch memory */
	complex_array scm;
	scm.resize(n + 2*m);
	
	int npeq= np+ 2*mp;
	nec_float fnop = nop;
	nec_float fnorm = 1.0/ fnop;
	int nrow= neq;
	
	if ( nop != 1)
	{
		for( ic = 0; ic < nrh; ic++ )
		{
			if ( (n != 0) && (m != 0) )
			{
				for( i = 0; i < neq; i++ )
					scm[i]= b[i+ic*neq];
			
				kk=2* mp;
				ia= np-1;
				ib= n-1;
				j= np-1;
			
				for( k = 0; k < nop; k++ )
				{
					if ( k != 0 )
					{
						for( i = 0; i < np; i++ )
						{
							ia++;
							j++;
							b[j+ic*neq]= scm[ia];
						}
				
						if ( k == (nop-1) )
							continue;
					} /* if ( k != 0 ) */
				
					for( i = 0; i < kk; i++ )
					{
						ib++;
						j++;
						b[j+ic*neq]= scm[ib];
					}
				} /* for( k = 0; k < nop; k++ ) */
			
			} /* if ( (n != 0) && (m != 0) ) */
		
			/* transform matrix eq. rhs vector according to symmetry modes */
			for( i = 0; i < npeq; i++ )
			{
				for( k = 0; k < nop; k++ )
				{
					ia= i+ k* npeq;
					scm[k]= b[ia+ic*neq];
				}
			
				sum= scm[0];
				for( k = 1; k < nop; k++ )
					sum += scm[k];
			
				b[i+ic*neq]= sum* fnorm;
			
				for( k = 1; k < nop; k++ )
				{
					ia= i+ k* npeq;
					sum= scm[0];
				
					for( j = 1; j < nop; j++ )
						sum += scm[j]* conj( symmetry_array[k+j*nop]);
				
					b[ia+ic*neq]= sum* fnorm;
				}
			} /* for( i = 0; i < npeq; i++ ) */
		
		} /* for( ic = 0; ic < nrh; ic++ ) */
	
	} /* if ( nop != 1) */
	
	/* solve each mode equation */
	for( kk = 0; kk < nop; kk++ )
	{
		ia= kk* npeq;
		ib= ia;
	
		for( ic = 0; ic < nrh; ic++ )
		{
			complex_array a_sub = a.sub_array(ib);
			complex_array b_sub = b.sub_array(ia+ic*neq);
			int_array ip_sub = ip.sub_array(ia);
			solve( npeq, a_sub, ip_sub, b_sub, nrow );
		}
	
	} /* for( kk = 0; kk < nop; kk++ ) */
	
	if ( nop == 1)
	{
		return;
	}
	
	/* inverse transform the mode solutions */
	for( ic = 0; ic < nrh; ic++ )
	{
		for( i = 0; i < npeq; i++ )
		{
			for( k = 0; k < nop; k++ )
			{
				ia= i+ k* npeq;
				scm[k]= b[ia+ic*neq];
			}
		
			sum= scm[0];
			for( k = 1; k < nop; k++ )
				sum += scm[k];
		
			b[i+ic*neq]= sum;
			for( k = 1; k < nop; k++ )
			{
				ia= i+ k* npeq;
				sum= scm[0];
			
				for( j = 1; j < nop; j++ )
					sum += scm[j]* symmetry_array[k+j*nop];
			
				b[ia+ic*neq]= sum;
			}
		
		} /* for( i = 0; i < npeq; i++ ) */
	
		if ( (n == 0) || (m == 0) )
			continue;
	
		for( i = 0; i < neq; i++ )
			scm[i]= b[i+ic*neq];
	
		kk=2* mp;
		ia= np-1;
		ib= n-1;
		j= np-1;
	
		for( k = 0; k < nop; k++ )
		{
			if ( k != 0 )
			{
				for( i = 0; i < np; i++ )
				{
					ia++;
					j++;
					b[ia+ic*neq]= scm[j];
				}
			
				if ( k == nop)
					continue;
			
			} /* if ( k != 0 ) */
		
			for( i = 0; i < kk; i++ )
			{
				ib++;
				j++;
				b[ib+ic*neq]= scm[j];
			}
		} /* for( k = 0; k < nop; k++ ) */
	
	} /* for( ic = 0; ic < nrh; ic++ ) */
}

/*-----------------------------------------------------------------------*/
/* test for convergence in numerical integration */
void test(
	nec_float f1r, nec_float f2r, nec_float *tr,
 	nec_float f1i, nec_float f2i, nec_float *ti,
	nec_float dmin )
{
	static nec_float _min_val =  1.0e-37;
	
/*
{
  double den;

  den= fabs( f2r);
  *tr= fabs( f2i);

  if( den < *tr)
    den= *tr;
  if( den < dmin)
    den= dmin;

  if( den < 1.0e-37)
  {
    *tr=0.;
    *ti=0.;
    return;
  }

  *tr= fabs(( f1r- f2r)/ den);
  *ti= fabs(( f1i- f2i)/ den);
}

*/	
	nec_float den = fabs( f2r);
	nec_float temp_tr = fabs( f2i);
	
	if( den < temp_tr)
		den = temp_tr;
	if( den < dmin)
		den = dmin;
	
	if( den < _min_val)
	{
		*tr = 0.0;
		*ti = 0.0;
		return;
	}
	
	*tr= fabs((f1r - f2r)/ den);
	*ti= fabs((f1i - f2i)/ den); 
}

/*
	Simpler test for convergence in numerical integration.
	This tests only one number. It is a special case of the
	test() function above.
*/
nec_float test_simple( nec_float f1r, nec_float f2r, nec_float dmin )
{
	static nec_float _min_val =  1.0e-37;
	
	nec_float den = fabs(f2r);
	
	if( den < dmin)
		den = dmin;
	if (den < _min_val)
	{
		return 0.0;
	}
	
	return fabs((f1r - f2r) / den);	
}