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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
/*================================================================
* mex_projection_QR.c = used by quickNcutHardBiais.m in eigensolver.
* mex_projection_QR(X,P,Ubar,R) == (eye(length(P))-Ubar*(R'*R)^{-1}*Ubar') * P * (eye(length(P))-Ubar*(R'*R)^{-1}*Ubar') * X ;
* (R'*R)^{-1} is solved by solving a triangular system
* P, Ubar, R are sparse, but X is full
* R is upper triangular
% test sequence:
*=================================================================*/
# include "math.h"
# include "mex.h"
# include "a_times_b_cmplx.c"
/*# include "a_times_b.c"*/
void mexFunction(
int nargout,
mxArray *out[],
int nargin,
const mxArray *in[]
)
{
int *ir[4], *jc[4], m[4], n[4];
double *y, *y1,*y2,*y3,*y4,*y5,*y6, *pr[4];
double *y2bis, *y5bis;
int k;
for (k=0; k<4; k++) {
m[k] = mxGetM(in[k]);
n[k] = mxGetN(in[k]);
pr[k] = mxGetPr(in[k]);
if (k>0) {
ir[k] = mxGetIr(in[k]);
jc[k] = mxGetJc(in[k]);
}
}
out[0] = mxCreateDoubleMatrix(m[1],1,mxREAL);
y = mxGetPr(out[0]);
y1 = mxCalloc(n[2], sizeof(double));
y2 = mxCalloc(m[3], sizeof(double));
y2bis = mxCalloc(m[3], sizeof(double));
y3 = mxCalloc(m[1], sizeof(double));
y4 = mxCalloc(m[1], sizeof(double));
y5 = mxCalloc(n[2], sizeof(double));
y5bis = mxCalloc(n[2], sizeof(double));
y6 = mxCalloc(n[2], sizeof(double));
CSR_VecMult_CAB_double(m[2],n[2],pr[2],ir[2],jc[2],pr[0],y1);
CSR_VecTriangSlvLD_CAB_double(m[3],pr[3],ir[3],jc[3],y1,y2bis);
CSC_VecTriangSlvUD_CAB_double(m[3],pr[3],ir[3],jc[3],y2bis,y2);
for(k=0;k<m[1];k++) {
y3[k]=pr[0][k];
}
CSC_VecMult_CaABC_double(m[2],n[2],-1,pr[2],ir[2],jc[2],y2,y3);
CSR_VecMult_CAB_double(m[1],n[1],pr[1],ir[1],jc[1],y3,y4);
CSR_VecMult_CAB_double(m[2],n[2],pr[2],ir[2],jc[2],y4,y5);
CSR_VecTriangSlvLD_CAB_double(m[3],pr[3],ir[3],jc[3],y5,y5bis);
CSC_VecTriangSlvUD_CAB_double(m[3],pr[3],ir[3],jc[3],y5bis,y6);
for(k=0;k<m[1];k++) {
y[k]=y4[k];
}
CSC_VecMult_CaABC_double(m[2],n[2],-1,pr[2],ir[2],jc[2],y6,y);
mxFree(y1);
mxFree(y2);
mxFree(y2bis);
mxFree(y3);
mxFree(y4);
mxFree(y5);
mxFree(y5bis);
mxFree(y6);
}
|
