本文整理汇总了C++中c_div函数的典型用法代码示例。如果您正苦于以下问题:C++ c_div函数的具体用法?C++ c_div怎么用?C++ c_div使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了c_div函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: small_conducting_Mie
void small_conducting_Mie(double x,struct c_complex m,double*mu,
long nangles,struct c_complex*s1,
struct c_complex*s2,double*qext,double*qsca,
double*qback,double*g)
/*:31*/
#line 448 "./mie.w"
{
struct c_complex ahat1,ahat2,bhat1,bhat2;
struct c_complex ss1;
double x2,x3,x4,muj,angle;
long j;
if((s1==NULL)||(s2==NULL))nangles= 0;
m.re+= 0.0;
x2= x*x;
x3= x2*x;
x4= x2*x2;
ahat1= c_div(c_set(0.0,2.0/3.0*(1.0-0.2*x2)),c_set(1.0-0.5*x2,2.0/3.0*x3));
bhat1= c_div(c_set(0.0,(x2-10.0)/30.0),c_set(1+0.5*x2,-x3/3.0));
ahat2= c_set(0.0,x2/30.);
bhat2= c_set(0.0,-x2/45.);
*qsca= 6.0*x4*(c_norm(ahat1)+c_norm(bhat1)+
(5.0/3.0)*(c_norm(ahat2)+c_norm(bhat2)));
*qext= *qsca;
*g= 6.0*x4*(ahat1.im*(ahat2.im+bhat1.im)+
bhat2.im*(5.0/9.0*ahat2.im+bhat1.im)+
ahat1.re*bhat1.re)/(*qsca);
ss1.re= 1.5*x2*(ahat1.re-bhat1.re);
ss1.im= 1.5*x2*(ahat1.im-bhat1.im-(5.0/3.0)*(ahat2.im+bhat2.im));
*qback= 4*c_norm(ss1);
x3*= 1.5;
ahat1.re*= x3;
ahat1.im*= x3;
bhat1.re*= x3;
bhat1.im*= x3;
ahat2.im*= x3*(5.0/3.0);
bhat2.im*= x3*(5.0/3.0);
for(j= 0;j<nangles;j++){
muj= mu[j];
angle= 2.0*muj*muj-1.0;
s1[j].re= ahat1.re+(bhat1.re)*muj;
s1[j].im= ahat1.im+(bhat1.im+ahat2.im)*muj+bhat2.im*angle;;
s2[j].re= bhat1.re+(ahat1.re)*muj;
s2[j].im= bhat1.im+(ahat1.im+bhat2.im)*muj+ahat2.im*angle;
}
}
开发者ID:JiapengHuang,项目名称:SPP,代码行数:54,代码来源:mie.c
示例2: Lentz_Dn
struct c_complex Lentz_Dn(struct c_complex z,long n)
/*:10*/
#line 126 "./mie.w"
{
struct c_complex alpha_j1,alpha_j2,zinv,aj;
struct c_complex alpha,result,ratio,runratio;
/*12:*/
#line 156 "./mie.w"
zinv= c_sdiv(2.0,z);
alpha= c_smul(n+0.5,zinv);
aj= c_smul(-n-1.5,zinv);
alpha_j1= c_add(aj,c_inv(alpha));
alpha_j2= aj;
ratio= c_div(alpha_j1,alpha_j2);
runratio= c_mul(alpha,ratio);
/*:12*/
#line 131 "./mie.w"
do
/*13:*/
#line 179 "./mie.w"
{
aj.re= zinv.re-aj.re;
aj.im= zinv.im-aj.im;
alpha_j1= c_add(c_inv(alpha_j1),aj);
alpha_j2= c_add(c_inv(alpha_j2),aj);
ratio= c_div(alpha_j1,alpha_j2);
zinv.re*= -1;
zinv.im*= -1;
runratio= c_mul(ratio,runratio);
}
/*:13*/
#line 134 "./mie.w"
while(fabs(c_abs(ratio)-1.0)> 1e-12);
result= c_add(c_sdiv((double)-n,z),runratio);
return result;
}
开发者ID:JiapengHuang,项目名称:SPP,代码行数:50,代码来源:mie.c
示例3: sur_print
static struct polynom *pol_div(struct polynom *p,struct polynom *q,struct polynom *p1,struct polynom *p2)
{
int i,j;
struct complex z,z1;
z.x=0; z.y=0; // RS 7.2.2013
p->n=p1->n-p2->n;
if (p->n<0)
{
sur_print("\nDegree of dividend < Degree of divisor");
WAIT; return(p);
}
for (i=0; i<=p->n; ++i)
p->a[i].x=p->a[i].y=0.0;
for (i=0; i<=p1->n; ++i)
{ q->a[i].x=p1->a[i].x; q->a[i].y=p1->a[i].y; }
for (i=p1->n; i>=p2->n; --i)
{
if (c_zero(&(q->a[i]))) continue;
c_div(&z,&(q->a[i]),&(p2->a[p2->n]));
p->a[i-p2->n].x=z.x; p->a[i-p2->n].y=z.y;
q->a[i].x=q->a[i].y=0.0;
for (j=i-1; j>=i-p2->n; --j)
c_sub(&(q->a[j]),&(q->a[j]),
c_mult(&z1,&z,&(p2->a[p2->n-i+j])));
}
i=p2->n-1;
while (c_zero(&(q->a[i])) && i>0) --i;
q->n=i;
return(p);
}
开发者ID:rforge,项目名称:muste,代码行数:32,代码来源:pol.c
示例4: expand_pole_product
/**
* \brief Expand pole product
* \param c resulting filter coefficients
* \param poles pole locations
* \param K number of poles
* \ingroup vyv_gaussian
*
* This routine expands the product to obtain the filter coefficients:
* \f[ \prod_{k=0}^{K-1}\frac{\mathrm{poles}[k]-1}{\mathrm{poles}[k]-z^{-1}}
= \frac{c[0]}{1+\sum_{k=1}^K c[k] z^{-k}}. \f]
*/
static void expand_pole_product(double *c, const complex4c *poles, int K)
{
complex4c denom[VYV_MAX_K + 1];
int k, j;
assert(K <= VYV_MAX_K);
denom[0] = poles[0];
denom[1] = make_complex(-1, 0);
for (k = 1; k < K; ++k)
{
denom[k + 1] = c_neg(denom[k]);
for (j = k; j > 0; --j)
denom[j] = c_sub(c_mul(denom[j], poles[k]), denom[j - 1]);
denom[0] = c_mul(denom[0], poles[k]);
}
for (k = 1; k <= K; ++k)
c[k] = c_div(denom[k], denom[0]).real;
for (c[0] = 1, k = 1; k <= K; ++k)
c[0] += c[k];
return;
}
开发者ID:ArtisticCoding,项目名称:OpenCP,代码行数:38,代码来源:gaussian_conv_vyv.cpp
示例5: dd
void dd (header *hd)
{ header *st=hd,*hd1,*result;
int c1,c2,i,j,r;
double *m1,*m2,*mr;
complex *mc1,*mc2,*mcr,hc1,hc2;
interval *mi1,*mi2,*mir,hi1,hi2;
hd1=next_param(st);
equal_params_2(&hd,&hd1); if (error) return;
getmatrix(hd,&r,&c1,&m1); if (r!=1) wrong_arg();
getmatrix(hd1,&r,&c2,&m2); if (r!=1) wrong_arg();
if (c1!=c2) wrong_arg();
if (iscomplex(hd)) /* complex values */
{ mc1=(complex *)m1; mc2=(complex *)m2;
result=new_cmatrix(1,c1,""); if (error) return;
mcr=(complex *)matrixof(result);
memmove((char *)mcr,(char *)mc2,c1*sizeof(complex));
for (i=1; i<c1; i++)
{ for (j=c1-1; j>=i; j--)
{ if (mc1[j][0]==mc1[j-i][0] &&
mc1[j][1]==mc1[j-i][1]) wrong_arg();
c_sub(mcr[j],mcr[j-1],hc1);
c_sub(mc1[j],mc1[j-i],hc2);
c_div(hc1,hc2,mcr[j]);
}
}
}
else if (isinterval(hd)) /* complex values */
{ mi1=(complex *)m1; mi2=(complex *)m2;
result=new_imatrix(1,c1,""); if (error) return;
mir=(interval *)matrixof(result);
memmove((char *)mir,(char *)mi2,c1*sizeof(interval));
for (i=1; i<c1; i++)
{ for (j=c1-1; j>=i; j--)
{ i_sub(mir[j],mir[j-1],hi1);
if (hi1[0]<=0 && hi1[1]>=0)
{ output("Interval points coincide\n");
error=1; return;
}
i_sub(mi1[j],mi1[j-i],hi2);
i_div(hi1,hi2,mir[j]);
}
}
}
else if (isreal(hd))
{ result=new_matrix(1,c1,""); if (error) return;
mr=matrixof(result);
memmove((char *)mr,(char *)m2,c1*sizeof(double));
for (i=1; i<c1; i++)
{ for (j=c1-1; j>=i; j--)
{ if (m1[j]==m1[j-i]) wrong_arg();
mr[j]=(mr[j]-mr[j-1])/(m1[j]-m1[j-i]);
}
}
}
else wrong_arg();
moveresult(st,result);
}
开发者ID:OS2World,项目名称:APP-MATH-Euler,代码行数:57,代码来源:polynom.cpp
示例6: variance
/**
* \brief Compute the variance of the impulse response
* \param poles0 unscaled pole locations
* \param q rescaling parameter
* \param K number of poles
* \return variance achieved by poles = poles0^(1/q)
* \ingroup vyv_gaussian
*/
static double variance(const complex4c *poles0, int K, double q)
{
complex4c sum = { 0, 0 };
int k;
for (k = 0; k < K; ++k)
{
complex4c z = c_real_pow(poles0[k], 1 / q), denom = z;
denom.real -= 1;
/* Compute sum += z / (z - 1)^2. */
sum = c_add(sum, c_div(z, c_mul(denom, denom)));
}
return 2 * sum.real;
}
开发者ID:ArtisticCoding,项目名称:OpenCP,代码行数:23,代码来源:gaussian_conv_vyv.cpp
示例7: dq_variance
/**
* \brief Derivative of variance with respect to q
* \param poles0 unscaled pole locations
* \param q rescaling parameter
* \param K number of poles
* \return derivative of variance with respect to q
* \ingroup vyv_gaussian
*
* This function is used by compute_q() in solving for q.
*/
static double dq_variance(const complex4c *poles0, int K, double q)
{
complex4c sum = { 0, 0 };
int k;
for (k = 0; k < K; ++k)
{
complex4c z = c_real_pow(poles0[k], 1 / q), w = z, denom = z;
w.real += 1;
denom.real -= 1;
/* Compute sum += z log(z) (z + 1) / (z - 1)^3 */
sum = c_add(sum, c_div(c_mul(c_mul(z, c_log(z)), w),
c_real_pow(denom, 3)));
}
return (2 / q) * sum.real;
}
开发者ID:ArtisticCoding,项目名称:OpenCP,代码行数:27,代码来源:gaussian_conv_vyv.cpp
示例8: cpivotL
//.........这里部分代码省略.........
int *lsub, *xlsub;
complex *lusup;
int *xlusup;
flops_t *ops = stat->ops;
/* Initialize pointers */
lsub = Glu->lsub;
xlsub = Glu->xlsub;
lusup = Glu->lusup;
xlusup = Glu->xlusup;
fsupc = (Glu->xsup)[(Glu->supno)[jcol]];
nsupc = jcol - fsupc; /* excluding jcol; nsupc >= 0 */
lptr = xlsub[fsupc];
nsupr = xlsub[fsupc+1] - lptr;
lu_sup_ptr = &lusup[xlusup[fsupc]]; /* start of the current supernode */
lu_col_ptr = &lusup[xlusup[jcol]]; /* start of jcol in the supernode */
lsub_ptr = &lsub[lptr]; /* start of row indices of the supernode */
#ifdef DEBUG
if ( jcol == MIN_COL ) {
printf(_T("Before cdiv: col %d\n"), jcol);
for (k = nsupc; k < nsupr; k++)
printf(_T(" lu[%d] %f\n"), lsub_ptr[k], lu_col_ptr[k]);
}
#endif
/* Determine the largest abs numerical value for partial pivoting;
Also search for user-specified pivot, and diagonal element. */
if ( *usepr ) *pivrow = iperm_r[jcol];
diagind = iperm_c[jcol];
pivmax = 0.0;
pivptr = nsupc;
diag = EMPTY;
old_pivptr = nsupc;
for (isub = nsupc; isub < nsupr; ++isub) {
rtemp = c_abs1 (&lu_col_ptr[isub]);
if ( rtemp > pivmax ) {
pivmax = rtemp;
pivptr = isub;
}
if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
if ( lsub_ptr[isub] == diagind ) diag = isub;
}
/* Test for singularity */
if ( pivmax == 0.0 ) {
*pivrow = lsub_ptr[pivptr];
perm_r[*pivrow] = jcol;
*usepr = 0;
return (jcol+1);
}
thresh = u * pivmax;
/* Choose appropriate pivotal element by our policy. */
if ( *usepr ) {
rtemp = c_abs1 (&lu_col_ptr[old_pivptr]);
if ( rtemp != 0.0 && rtemp >= thresh )
pivptr = old_pivptr;
else
*usepr = 0;
}
if ( *usepr == 0 ) {
/* Use diagonal pivot? */
if ( diag >= 0 ) { /* diagonal exists */
rtemp = c_abs1 (&lu_col_ptr[diag]);
if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = diag;
}
*pivrow = lsub_ptr[pivptr];
}
/* Record pivot row */
perm_r[*pivrow] = jcol;
/* Interchange row subscripts */
if ( pivptr != nsupc ) {
itemp = lsub_ptr[pivptr];
lsub_ptr[pivptr] = lsub_ptr[nsupc];
lsub_ptr[nsupc] = itemp;
/* Interchange numerical values as well, for the whole snode, such
* that L is indexed the same way as A.
*/
for (icol = 0; icol <= nsupc; icol++) {
itemp = pivptr + icol * nsupr;
temp = lu_sup_ptr[itemp];
lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];
lu_sup_ptr[nsupc + icol*nsupr] = temp;
}
} /* if */
/* cdiv operation */
ops[FACT] += 10 * (nsupr - nsupc);
c_div(&temp, &one, &lu_col_ptr[nsupc]);
for (k = nsupc+1; k < nsupr; k++)
cc_mult(&lu_col_ptr[k], &lu_col_ptr[k], &temp);
return 0;
}
开发者ID:artemeliy,项目名称:inf4715,代码行数:101,代码来源:cpivotL.c
示例9: lsame_
/* Subroutine */ int ctrti2_(char *uplo, char *diag, integer *n, complex *a,
integer *lda, integer *info)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
complex q__1;
/* Local variables */
integer j;
complex ajj;
logical upper;
logical nounit;
/* -- LAPACK routine (version 3.2) -- */
/* November 2006 */
/* Purpose */
/* ======= */
/* CTRTI2 computes the inverse of a complex upper or lower triangular */
/* matrix. */
/* This is the Level 2 BLAS version of the algorithm. */
/* Arguments */
/* ========= */
/* UPLO (input) CHARACTER*1 */
/* Specifies whether the matrix A is upper or lower triangular. */
/* = 'U': Upper triangular */
/* = 'L': Lower triangular */
/* DIAG (input) CHARACTER*1 */
/* Specifies whether or not the matrix A is unit triangular. */
/* = 'N': Non-unit triangular */
/* = 'U': Unit triangular */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* A (input/output) COMPLEX array, dimension (LDA,N) */
/* On entry, the triangular matrix A. If UPLO = 'U', the */
/* leading n by n upper triangular part of the array A contains */
/* the upper triangular matrix, and the strictly lower */
/* triangular part of A is not referenced. If UPLO = 'L', the */
/* leading n by n lower triangular part of the array A contains */
/* the lower triangular matrix, and the strictly upper */
/* triangular part of A is not referenced. If DIAG = 'U', the */
/* diagonal elements of A are also not referenced and are */
/* assumed to be 1. */
/* On exit, the (triangular) inverse of the original matrix, in */
/* the same storage format. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,N). */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -k, the k-th argument had an illegal value */
/* ===================================================================== */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
*info = 0;
upper = lsame_(uplo, "U");
nounit = lsame_(diag, "N");
if (! upper && ! lsame_(uplo, "L")) {
*info = -1;
} else if (! nounit && ! lsame_(diag, "U")) {
*info = -2;
} else if (*n < 0) {
*info = -3;
} else if (*lda < max(1,*n)) {
*info = -5;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("CTRTI2", &i__1);
return 0;
}
if (upper) {
/* Compute inverse of upper triangular matrix. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (nounit) {
i__2 = j + j * a_dim1;
c_div(&q__1, &c_b1, &a[j + j * a_dim1]);
a[i__2].r = q__1.r, a[i__2].i = q__1.i;
i__2 = j + j * a_dim1;
//.........这里部分代码省略.........
开发者ID:juanjosegarciaripoll,项目名称:cblapack,代码行数:101,代码来源:ctrti2.c
示例10: twonorm
/* Subroutine */ int claic1_(integer *job, integer *j, complex *x, real *sest,
complex *w, complex *gamma, real *sestpr, complex *s, complex *c__)
{
/* -- LAPACK auxiliary routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
October 31, 1992
Purpose
=======
CLAIC1 applies one step of incremental condition estimation in
its simplest version:
Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
lower triangular matrix L, such that
twonorm(L*x) = sest
Then CLAIC1 computes sestpr, s, c such that
the vector
[ s*x ]
xhat = [ c ]
is an approximate singular vector of
[ L 0 ]
Lhat = [ w' gamma ]
in the sense that
twonorm(Lhat*xhat) = sestpr.
Depending on JOB, an estimate for the largest or smallest singular
value is computed.
Note that [s c]' and sestpr**2 is an eigenpair of the system
diag(sest*sest, 0) + [alpha gamma] * [ conjg(alpha) ]
[ conjg(gamma) ]
where alpha = conjg(x)'*w.
Arguments
=========
JOB (input) INTEGER
= 1: an estimate for the largest singular value is computed.
= 2: an estimate for the smallest singular value is computed.
J (input) INTEGER
Length of X and W
X (input) COMPLEX array, dimension (J)
The j-vector x.
SEST (input) REAL
Estimated singular value of j by j matrix L
W (input) COMPLEX array, dimension (J)
The j-vector w.
GAMMA (input) COMPLEX
The diagonal element gamma.
SESTPR (output) REAL
Estimated singular value of (j+1) by (j+1) matrix Lhat.
S (output) COMPLEX
Sine needed in forming xhat.
C (output) COMPLEX
Cosine needed in forming xhat.
=====================================================================
Parameter adjustments */
/* Table of constant values */
static integer c__1 = 1;
/* System generated locals */
real r__1, r__2;
complex q__1, q__2, q__3, q__4, q__5, q__6;
/* Builtin functions */
double c_abs(complex *);
void r_cnjg(complex *, complex *), c_sqrt(complex *, complex *);
double sqrt(doublereal);
void c_div(complex *, complex *, complex *);
/* Local variables */
static complex sine;
static real test, zeta1, zeta2, b, t;
static complex alpha;
extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer
*, complex *, integer *);
static real norma, s1, s2, absgam, absalp;
extern doublereal slamch_(char *);
static complex cosine;
static real absest, scl, eps, tmp;
--w;
--x;
/* Function Body */
//.........这里部分代码省略.........
开发者ID:EugeneGalipchak,项目名称:antelope_contrib,代码行数:101,代码来源:claic1.c
示例11: r_imag
/* Subroutine */ int cgttrf_(integer *n, complex *dl, complex *d__, complex *
du, complex *du2, integer *ipiv, integer *info)
{
/* System generated locals */
integer i__1, i__2, i__3, i__4;
real r__1, r__2, r__3, r__4;
complex q__1, q__2;
/* Builtin functions */
double r_imag(complex *);
void c_div(complex *, complex *, complex *);
/* Local variables */
integer i__;
complex fact, temp;
extern /* Subroutine */ int xerbla_(char *, integer *);
/* -- LAPACK routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CGTTRF computes an LU factorization of a complex tridiagonal matrix A */
/* using elimination with partial pivoting and row interchanges. */
/* The factorization has the form */
/* A = L * U */
/* where L is a product of permutation and unit lower bidiagonal */
/* matrices and U is upper triangular with nonzeros in only the main */
/* diagonal and first two superdiagonals. */
/* Arguments */
/* ========= */
/* N (input) INTEGER */
/* The order of the matrix A. */
/* DL (input/output) COMPLEX array, dimension (N-1) */
/* On entry, DL must contain the (n-1) sub-diagonal elements of */
/* A. */
/* On exit, DL is overwritten by the (n-1) multipliers that */
/* define the matrix L from the LU factorization of A. */
/* D (input/output) COMPLEX array, dimension (N) */
/* On entry, D must contain the diagonal elements of A. */
/* On exit, D is overwritten by the n diagonal elements of the */
/* upper triangular matrix U from the LU factorization of A. */
/* DU (input/output) COMPLEX array, dimension (N-1) */
/* On entry, DU must contain the (n-1) super-diagonal elements */
/* of A. */
/* On exit, DU is overwritten by the (n-1) elements of the first */
/* super-diagonal of U. */
/* DU2 (output) COMPLEX array, dimension (N-2) */
/* On exit, DU2 is overwritten by the (n-2) elements of the */
/* second super-diagonal of U. */
/* IPIV (output) INTEGER array, dimension (N) */
/* The pivot indices; for 1 <= i <= n, row i of the matrix was */
/* interchanged with row IPIV(i). IPIV(i) will always be either */
/* i or i+1; IPIV(i) = i indicates a row interchange was not */
/* required. */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -k, the k-th argument had an illegal value */
/* > 0: if INFO = k, U(k,k) is exactly zero. The factorization */
/* has been completed, but the factor U is exactly */
/* singular, and division by zero will occur if it is used */
/* to solve a system of equations. */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Statement Functions .. */
/* .. */
/* .. Statement Function definitions .. */
/* .. */
/* .. Executable Statements .. */
/* Parameter adjustments */
//.........这里部分代码省略.........
开发者ID:Gaylou,项目名称:CMVS-PMVS,代码行数:101,代码来源:cgttrf.c
示例12: r_imag
/*< subroutine csvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) >*/
/* Subroutine */ int csvdc_(complex *x, integer *ldx, integer *n, integer *p,
complex *s, complex *e, complex *u, integer *ldu, complex *v, integer
*ldv, complex *work, integer *job, integer *info)
{
/* System generated locals */
integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2,
i__3, i__4;
real r__1, r__2, r__3, r__4;
complex q__1, q__2, q__3;
/* Builtin functions */
double r_imag(complex *), c_abs(complex *);
void c_div(complex *, complex *, complex *), r_cnjg(complex *, complex *);
double sqrt(doublereal);
/* Local variables */
real b, c__, f, g;
integer i__, j, k, l=0, m;
complex r__, t;
real t1, el;
integer kk;
real cs;
integer ll, mm, ls=0;
real sl;
integer lu;
real sm, sn;
integer lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
real emm1, smm1;
integer kase, jobu, iter;
real test;
integer nctp1, nrtp1;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *);
real scale;
extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer
*, complex *, integer *);
real shift;
extern /* Subroutine */ int cswap_(integer *, complex *, integer *,
complex *, integer *);
integer maxit;
extern /* Subroutine */ int caxpy_(integer *, complex *, complex *,
integer *, complex *, integer *), csrot_(integer *, complex *,
integer *, complex *, integer *, real *, real *);
logical wantu, wantv;
extern /* Subroutine */ int srotg_(real *, real *, real *, real *);
real ztest;
extern doublereal scnrm2_(integer *, complex *, integer *);
/*< integer ldx,n,p,ldu,ldv,job,info >*/
/*< complex x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1) >*/
/* csvdc is a subroutine to reduce a complex nxp matrix x by */
/* unitary transformations u and v to diagonal form. the */
/* diagonal elements s(i) are the singular values of x. the */
/* columns of u are the corresponding left singular vectors, */
/* and the columns of v the right singular vectors. */
/* on entry */
/* x complex(ldx,p), where ldx.ge.n. */
/* x contains the matrix whose singular value */
/* decomposition is to be computed. x is */
/* destroyed by csvdc. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix x. */
/* p integer. */
/* p is the number of columns of the matrix x. */
/* ldu integer. */
/* ldu is the leading dimension of the array u */
/* (see below). */
/* ldv integer. */
/* ldv is the leading dimension of the array v */
/* (see below). */
/* work complex(n). */
/* work is a scratch array. */
/* job integer. */
/* job controls the computation of the singular */
/* vectors. it has the decimal expansion ab */
/* with the following meaning */
/* a.eq.0 do not compute the left singular */
/* vectors. */
/* a.eq.1 return the n left singular vectors */
/* in u. */
/* a.ge.2 returns the first min(n,p) */
/* left singular vectors in u. */
/* b.eq.0 do not compute the right singular */
/* vectors. */
/* b.eq.1 return the right singular vectors */
//.........这里部分代码省略.........
开发者ID:151706061,项目名称:ITK,代码行数:101,代码来源:csvdc.c
示例13: sqrt
/** CHETRF_ROOK_REC2 computes a partial factorization of a complex Hermitian indefinite matrix using the boun ded Bunch-Kaufman ("rook") diagonal pivoting method
*
* This routine is a minor modification of LAPACK's clahef_rook.
* It serves as an unblocked kernel in the recursive algorithms.
* The blocked BLAS Level 3 updates were removed and moved to the
* recursive algorithm.
* */
/* Subroutine */ void RELAPACK_chetrf_rook_rec2(char *uplo, int *n,
int *nb, int *kb, complex *a, int *lda, int *ipiv,
complex *w, int *ldw, int *info, ftnlen uplo_len)
{
/* System generated locals */
int a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3, i__4;
float r__1, r__2;
complex q__1, q__2, q__3, q__4, q__5;
/* Builtin functions */
double sqrt(double), r_imag(complex *);
void r_cnjg(complex *, complex *), c_div(complex *, complex *, complex *);
/* Local variables */
static int j, k, p;
static float t, r1;
static complex d11, d21, d22;
static int ii, jj, kk, kp, kw, jp1, jp2, kkw;
static logical done;
static int imax, jmax;
static float alpha;
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern /* Subroutine */ int cgemv_(char *, int *, int *, complex *
, complex *, int *, complex *, int *, complex *, complex *
, int *, ftnlen);
static float sfmin;
extern /* Subroutine */ int ccopy_(int *, complex *, int *,
complex *, int *);
static int itemp;
extern /* Subroutine */ int cswap_(int *, complex *, int *,
complex *, int *);
static int kstep;
static float stemp, absakk;
extern /* Subroutine */ int clacgv_(int *, complex *, int *);
extern int icamax_(int *, complex *, int *);
extern double slamch_(char *, ftnlen);
extern /* Subroutine */ int csscal_(int *, float *, complex *, int
*);
static float colmax, rowmax;
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--ipiv;
w_dim1 = *ldw;
w_offset = 1 + w_dim1;
w -= w_offset;
/* Function Body */
*info = 0;
alpha = (sqrt(17.f) + 1.f) / 8.f;
sfmin = slamch_("S", (ftnlen)1);
if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
k = *n;
L10:
kw = *nb + k - *n;
if ((k <= *n - *nb + 1 && *nb < *n) || k < 1) {
goto L30;
}
kstep = 1;
p = k;
if (k > 1) {
i__1 = k - 1;
ccopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &w[kw * w_dim1 + 1], &
c__1);
}
i__1 = k + kw * w_dim1;
i__2 = k + k * a_dim1;
r__1 = a[i__2].r;
w[i__1].r = r__1, w[i__1].i = 0.f;
if (k < *n) {
i__1 = *n - k;
q__1.r = -1.f, q__1.i = -0.f;
cgemv_("No transpose", &k, &i__1, &q__1, &a[(k + 1) * a_dim1 + 1],
lda, &w[k + (kw + 1) * w_dim1], ldw, &c_b1, &w[kw *
w_dim1 + 1], &c__1, (ftnlen)12);
i__1 = k + kw * w_dim1;
i__2 = k + kw * w_dim1;
r__1 = w[i__2].r;
w[i__1].r = r__1, w[i__1].i = 0.f;
}
i__1 = k + kw * w_dim1;
absakk = (r__1 = w[i__1].r, dabs(r__1));
if (k > 1) {
i__1 = k - 1;
imax = icamax_(&i__1, &w[kw * w_dim1 + 1], &c__1);
i__1 = imax + kw * w_dim1;
colmax = (r__1 = w[i__1].r, dabs(r__1)) + (r__2 = r_imag(&w[imax
+ kw * w_dim1]), dabs(r__2));
} else {
colmax = 0.f;
}
//.........这里部分代码省略.........
开发者ID:ChinaQuants,项目名称:OpenBLAS,代码行数:101,代码来源:chetrf_rook_rec2.c
示例14: polroot
static int polroot(
struct polynom *p,
struct complex *pz, /* pointer to root */
struct complex *pz0 /* pointer to initial value */
)
{
int i;
int n_iter;
struct polynom d;
struct complex v,v0,vd,delta,zmin;
double y,ymin;
// extern struct polynom *pol_der();
// extern struct complex *pol_value();
delta.x=0; delta.y=0; // RS 7.2.2013
//printf("\np->n=%d",p->n); getch();
if (p->n==1)
{
c_div(pz,&(p->a[0]),&(p->a[1]));
pz->x=-(pz->x); pz->y=-(pz->y);
return(1);
}
//printf("\nr2");
pol_value(p,pz0,&v0);
//printf("\nr3");
if (c_zero(&v0))
{
pz->x=pz0->x;
pz->y=pz0->y;
return(1);
}
//printf("\nr4");
zmin.x=pz0->x; zmin.y=pz0->y;
ymin=v0.x*v0.x+v0.y*v0.y;
pol_der(&d,p);
//printf("\nr5");
n_iter=0;
while (1)
{
pol_value(&d,pz0,&vd);
c_div(&delta,&v0,&vd);
c_sub(pz,pz0,&delta);
pol_value(p,pz,&v);
v0.x=v.x; v0.y=v.y;
pz0->x=pz->x; pz0->y=pz->y;
++n_iter;
PR_UP;
sprintf(sbuf,"\npolroot: N=%d Re=%e Im=%e",n_iter,pz->x,pz->y); sur_print(sbuf); // RS CHA Rprintf
/* Rprintf("zero: %d\n",c_zero(pz)); getch(); */
/* if (c_zero(&pz)) break; */
if (c_zero(pz)) break;
y=v.x*v.x+v.y*v.y;
if (y<ymin)
{ ymin=y; zmin.x=pz->x; zmin.y=pz->y; }
if (fabs(delta.x)<roots_eps && fabs(delta.y)<roots_eps) break;
if (n_iter>roots_max_iter)
{
pz->x=zmin.x; pz->y=zmin.y;
break;
}
if (sur_kbhit())
{
pz->x=zmin.x; pz->y=zmin.y;
i=sur_getch(); if (i=='.') break;
}
}
return (n_iter);
}
开发者ID:rforge,项目名称:muste,代码行数:72,代码来源:pol.c
示例15: c_abs
/* Subroutine */ int clarge_(integer *n, complex *a, integer *lda, integer *
iseed, complex *work, integer *info)
{
/* System generated locals */
integer a_dim1, a_offset, i__1;
real r__1;
complex q__1;
/* Builtin functions */
double c_abs(complex *);
void c_div(complex *, complex *, complex *);
/* Local variables */
static integer i__;
extern /* Subroutine */ int cgerc_(integer *, integer *, complex *,
complex *, integer *, complex *, integer *, complex *, integer *),
cscal_(integer *, complex *, complex *, integer *), cgemv_(char *
, integer *, integer *, complex *, complex *, integer *, complex *
, integer *, complex *, complex *, integer *);
extern doublereal scnrm2_(integer *, complex *, integer *);
static complex wa, wb;
static real wn;
extern /* Subroutine */ int xerbla_(char *, integer *), clarnv_(
integer *, integer *, integer *, complex *);
static complex tau;
#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1
#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]
/* -- LAPACK auxiliary test routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
September 30, 1994
Purpose
=======
CLARGE pre- and post-multiplies a complex general n by n matrix A
with a random unitary matrix: A = U*D*U'.
Arguments
=========
N (input) INTEGER
The order of the matrix A. N >= 0.
A (input/output) COMPLEX array, dimension (LDA,N)
On entry, the original n by n matrix A.
On exit, A is overwritten by U*A*U' for some random
unitary matrix U.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= N.
ISEED (input/output) INTEGER array, dimension (4)
On entry, the seed of the random number generator; the array
elements must be between 0 and 4095, and ISEED(4) must be
odd.
On exit, the seed is updated.
WORK (workspace) COMPLEX array, dimension (2*N)
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
=====================================================================
Test the input arguments
Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
--iseed;
--work;
/* Function Body */
*info = 0;
if (*n < 0) {
*info = -1;
} else if (*lda < max(1,*n)) {
*info = -3;
}
if (*info < 0) {
i__1 = -(*info);
xerbla_("CLARGE", &i__1);
return 0;
}
/* pre- and post-multiply A by random unitary matrix */
for (i__ = *n; i__ >= 1; --i__) {
/* generate random reflection */
//.........这里部分代码省略.........
开发者ID:zangel,项目名称:uquad,代码行数:101,代码来源:clarge.c
示例16: norm
//.........这里部分代码省略.........
} else if (normdif != 0.f) {
nwise_err__ = slamch_("OVERFLOW");
} else {
nwise_err__ = 0.f;
}
i__4 = n;
for (i__ = 1; i__ <= i__4; ++i__) {
rinv[i__ - 1] = 0.f;
}
i__4 = n;
for (j = 1; j <= i__4; ++j) {
i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
i__3 = i__ + j * 6 - 7;
i__5 = j + k * 6 - 7;
q__2.r = a[i__3].r * invhilb[i__5].r - a[i__3].i *
invhilb[i__5].i, q__2.i = a[i__3].r * invhilb[
i__5].i + a[i__3].i * invhilb[i__5].r;
q__1.r = q__2.r, q__1.i = q__2.i;
rinv[i__ - 1] += (r__1 = q__1.r, dabs(r__1)) + (r__2 =
r_imag(&q__1), dabs(r__2));
}
}
rinorm = 0.f;
i__4 = n;
for (i__ = 1; i__ <= i__4; ++i__) {
sumri = 0.f;
i__2 = n;
for (j = 1; j <= i__2; ++j) {
i__3 = i__ + j * 6 - 7;
i__5 = j - 1;
q__3.r = rinv[i__5] * invhilb[i__3].r, q__3.i = rinv[i__5]
* invhilb[i__3].i;
c_div(&q__2, &q__3, &invhilb[i__ + k * 6 - 7]);
q__1.r = q__2.r, q__1.i = q__2.i;
sumri += (r__1 = q__1.r, dabs(r__1)) + (r__2 = r_imag(&
q__1), dabs(r__2));
}
rinorm = dmax(rinorm,sumri);
}
/* invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm */
/* by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) */
ccond = ((r__1 = a[0].r, dabs(r__1)) + (r__2 = r_imag(a), dabs(
r__2))) / rinorm;
/* Forward error bound tests */
nwise_bnd__ = errbnd_n__[k + nrhs - 1];
cwise_bnd__ = errbnd_c__[k + nrhs - 1];
nwise_rcond__ = errbnd_n__[k + (nrhs << 1) - 1];
cwise_rcond__ = errbnd_c__[k + (nrhs << 1) - 1];
/* write (*,*) 'nwise : ', n, k, ncond, nwise_rcond, */
/* $ condthresh, ncond.ge.condthresh */
/* write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh */
if (ncond >= condthresh) {
s_copy(nguar, "YES", (ftnlen)3, (ftnlen)3);
if (nwise_bnd__ > errthresh) {
tstrat[0] = 1 / (eps * 2.f);
} else {
if (nwise_bnd__ != 0.f) {
tstrat[0] = nwise_err__ / nwise_bnd__;
} else if (nwise_err__ != 0.f) {
tstrat[0] = 1 / (eps * 16.f);
} else {
tstrat[0] = 0.f;
}
if (tstrat[0] > 1.f) {
tstrat[0] = 1 / (eps * 4.f);
开发者ID:juanjosegarciaripoll,项目名称:cblapack,代码行数:67,代码来源:cebchvxx.c
示例17: c_div
/* Subroutine */ int csytrs_(char *uplo, integer *n, integer *nrhs, complex *
a, integer *lda, integer *ipiv, complex *b, integer *ldb, integer *
info)
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
complex q__1, q__2, q__3;
/* Builtin functions */
void c_div(complex *, complex *, complex *);
/* Local variables */
integer j, k;
complex ak, bk;
integer kp;
complex akm1, bkm1, akm1k;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *);
extern logical lsame_(char *, char *);
complex denom;
extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
, complex *, integer *, complex *, integer *, complex *, complex *
, integer *), cgeru_(integer *, integer *, complex *,
complex *, integer *, complex *, integer *, complex *, integer *),
cswap_(integer *, complex *, integer *, complex *, integer *);
logical upper;
extern /* Subroutine */ int xerbla_(char *, integer *);
/* -- LAPACK routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CSYTRS solves a system of linear equations A*X = B with a complex */
/* symmetric matrix A using the factorization A = U*D*U**T or */
/* A = L*D*L**T computed by CSYTRF. */
/* Arguments */
/* ========= */
/* UPLO (input) CHARACTER*1 */
/* Specifies whether the details of the factorization are stored */
/* as an upper or lower triangular matrix. */
/* = 'U': Upper triangular, form is A = U*D*U**T; */
/* = 'L': Lower triangular, form is A = L*D*L**T. */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* NRHS (input) INTEGER */
/* The number of right hand sides, i.e., the number of columns */
/* of the matrix B. NRHS >= 0. */
/* A (input) COMPLEX array, dimension (LDA,N) */
/* The block diagonal matrix D and the multipliers used to */
/* obtain the factor U or L as computed by CSYTRF. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,N). */
/* IPIV (input) INTEGER array, dimension (N) */
/* Details of the interchanges and the block structure of D */
/* as determined by CSYTRF. */
/* B (input/output) COMPLEX array, dimension (LDB,NRHS) */
/* On entry, the right hand side matrix B. */
/* On exit, the solution matrix X. */
/* LDB (input) INTEGER */
/* The leading dimension of the array B. LDB >= max(1,N). */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
//.........这里部分代码省略.........
开发者ID:3deggi,项目名称:levmar-ndk,代码行数:101,代码来源:csytrs.c
示例18: ilu_cpivotL
//.........这里部分代码省略.........
else {
thresh = u * pivmax;
/* Choose appropriate pivotal element by our policy. */
if ( *usepr ) {
switch (milu) {
case SMILU_1:
c_add(&temp, &lu_col_ptr[old_pivptr], &drop_sum);
rtemp = c_abs1(&temp);
break;
case SMILU_2:
case SMILU_3:
rtemp = c_abs1(&lu_col_ptr[old_pivptr]) + drop_sum.r;
break;
case SILU:
|
请发表评论