qr.c

最大似然估计算法

#include "timeseries.h"

void qr(double *A, int M, int N, double *Q, double *R) 
{

int j, k, l;
int ntau, info, lwork;

double *x, *tau, *work;

ntau = M < N ? M : N;
lwork = 2 * N;

x    = (double *) calloc((size_t)(M*N),  sizeof(double));
tau  = (double *) calloc((size_t) ntau,  sizeof(double));
work = (double *) calloc((size_t) lwork, sizeof(double));

dlacpy_("Full", &M, &N, A, &M, x, &M);

dgeqrf_(&M,&N,x,&M,tau,work,&lwork,&info);

for (j = 0; j < ntau; j++) {
	for (k = 0; k < N; k++) R[j + k * ntau] = 0.0;
}
for (j = 0; j < ntau; j++) {
	for (k = j; k < N; k++) R[j + k * ntau] = x[j + k * M];
}

dorgqr_(&M,&N,&N,x,&M,tau,work,&lwork,&info);

for (j = 0; j < M; j++) {
	for (k = 0; k < ntau; k++) Q[j + k * M] = x[j+ k * M];
}

free(x);
free(tau);
free(work);

}