1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
| void qrdcmp(float **A, int N, float *c, float *d, int *sing)
{
int i,j,k;
float scale,sigma,sum,tau;
*sing=0;
for (k=0;k<N-1;k++)
{
scale=0.0;
for (i=k;i<N;i++)
{
if (scale>=fabs(A[i][k]))
{
scale=scale;
}
else
{
scale=fabs(A[i][k]);
}
}
if (scale == 0.0)
{
*sing=1;
c[k]=d[k]=0.0;
}
else
{
for (i=k;i<N;i++) A[i][k] /= scale;
for (sum=0.0,i=k;i<=N;i++) sum += (A[i][k])*(A[i][k]); if (A[k][k]>=0)
{
sigma=sqrt(sum);
}
else
{
sigma=-sqrt(sum);
}
A[k][k] += sigma;
c[k]=sigma*A[k][k];
d[k] = -scale*sigma;
for (j=k+1;j<N;j++)
{
for (sum=0.0,i=k;i<=N;i++) sum += A[i][k]*A[i][j]; tau=sum/c[k];
for (i=k;i<N;i++) A[i][j] -= tau*A[i][k];
}
}
}
d[N]=A[N-1][N-1];
if (d[N-1] == 0.0) *sing=1;
} |
Partager