20static int SMatDestroy(
void*S){
25 DSDPFREE(&SS->sinv,&info);
31static int SMatGetSize(
void *S,
int *n){
37static int SMatView(
void* S){
40 info=Mat4View(SS->spsym); DSDPCHKERR(info);
44static int SMatLogDet(
void* S,
double *dd){
47 info=Mat4LogDet(SS->spsym,dd); DSDPCHKERR(info);
53static int SMatSetURMatP(
void*S,
double v[],
int nn,
int n){
56 double *rw1,*rw2,*xr=v;
57 rw1=SS->spsym->rw;rw2=rw1+n;
58 info=MatZeroEntries4(SS->spsym); DSDPCHKERR(info);
63 memcpy((
void*)rw1,(
void*)(xr),(row+1)*
sizeof(
double));
66 memcpy((
void*)(rw2),(
void*)(xr),(row+2)*
sizeof(
double));
70 for (j=row+2;j<n;j++){
76 info=MatSetColumn4(SS->spsym,rw1,row); DSDPCHKERR(info);
77 info=MatSetColumn4(SS->spsym,rw2,row+1); DSDPCHKERR(info);
80 for (row=2*(n/2);row<n;row++){
83 memcpy((
void*)(rw1),(
void*)(xr),(row+1)*
sizeof(
double));
85 for (j=row+1;j<n;j++){ rw1[j]=xr[row]; xr+=(j+2); }
86 info=MatSetColumn4(SS->spsym,rw1,row); DSDPCHKERR(info);
92static int SMatSetURMatU(
void*S,
double v[],
int nn,
int n){
95 double *rw1,*rw2,*xr=v;
96 rw1=SS->spsym->rw;rw2=rw1+n;
97 info=MatZeroEntries4(SS->spsym); DSDPCHKERR(info);
102 memcpy((
void*)rw1,(
void*)(xr),(row+1)*
sizeof(
double));
105 memcpy((
void*)(rw2),(
void*)(xr),(row+2)*
sizeof(
double));
108 for (j=row+2;j<n;j++){
114 info=MatSetColumn4(SS->spsym,rw1,row); DSDPCHKERR(info);
115 info=MatSetColumn4(SS->spsym,rw2,row+1); DSDPCHKERR(info);
118 for (row=2*(n/2);row<n;row++){
121 memcpy((
void*)(rw1),(
void*)(xr),(row+1)*
sizeof(
double));
123 for (j=row+1;j<n;j++){ rw1[j]=xr[row]; xr+=n; }
124 info=MatSetColumn4(SS->spsym,rw1,row); DSDPCHKERR(info);
129static int SMatSetURMat(
void*S,
double v[],
int nn,
int n){
133 info=SMatSetURMatP(S,v,nn,n);DSDPCHKERR(info);
134 }
else if (SS->UPLQ==
'U'){
135 info=SMatSetURMatU(S,v,nn,n);DSDPCHKERR(info);
140static int SMatSolve(
void *S,
int indx[],
int nind,
double b[],
double x[],
int n){
143 double alpha,*s1=SS->sinv,*s2;
145 if (SS->sinv && nind < n/4){
146 memset((
void*)x,0,n*
sizeof(
double));
147 for (i=0;i<nind;i++){
149 ione=1;nn=n;alpha=b[ii];s2=s1+n*ii;
150 daxpy(&nn,&alpha,s2,&ione,x,&ione);
153 memcpy(x,b,n*
sizeof(
double));
154 ChlSolve(SS->spsym, b, x);
159static int SMatCholeskySolveBackward(
void *S,
double b[],
double x[],
int n){
161 ChlSolveBackward(SS->spsym, b, x);
165static int SMatCholeskySolveForward(
void *S,
double b[],
double x[],
int n){
167 ChlSolveForward(SS->spsym, b, x);
171static int SMatFull(
void *S,
int *full){
176static int SMatCholeskyFactor(
void *S,
int *flag){
183 Cfact=(cfc_sta)ChlFact(SS->spsym,iw,rw,TRUE);
192static int SMatInverseAddP(
void *S,
double alpha,
double v[],
int nn,
int n){
195 double *x,*b,*ss=SS->sinv;
201 daxpy(&ii,&alpha,ss,&ione,v,&ione);
205 b=SS->spsym->rw;x=b+n;
207 memset((
void*)b,0,n*
sizeof(
double));
209 ChlSolve(SS->spsym, b, x);
219static int SMatInverseAddU(
void *S,
double alpha,
double v[],
int nn,
int n){
221 ffinteger n2=n*n,ione=1;
222 double *x,*b,*ss=SS->sinv;
226 daxpy(&n2,&alpha,ss,&ione,v,&ione);
228 b=SS->spsym->rw;x=b+n;
230 memset((
void*)b,0,n*
sizeof(
double));
232 ChlSolve(SS->spsym, b, x);
242static int SMatInverseAdd(
void *S,
double alpha,
double v[],
int nn,
int n){
246 info=SMatInverseAddP(S,alpha,v,nn,n);DSDPCHKERR(info);
247 }
else if (SS->UPLQ==
'U'){
248 info=SMatInverseAddU(S,alpha,v,nn,n);DSDPCHKERR(info);
253static int SMatCholeskyForwardMultiply(
void *S,
double b[],
double x[],
int n){
255 GetUhat(SS->spsym,b,x);
259static int SMatInvert(
void *S){
261 double *x,*b,*v=SS->sinv;
264 b=SS->spsym->rw;x=b+n;
268 memset((
void*)b,0,n*
sizeof(
double));
270 ChlSolve(SS->spsym, b, x);
271 memcpy((
void*)(v+i*n),(
void*)x,n*
sizeof(
double));
278static const char* tmatname=
"SPARSE PSD";
281 if (sops==NULL)
return 0;
283 sops->matcholesky=SMatCholeskyFactor;
284 sops->matsolveforward=SMatCholeskySolveForward;
285 sops->matsolvebackward=SMatCholeskySolveBackward;
286 sops->matinversemultiply=SMatSolve;
287 sops->matinvert=SMatInvert;
288 sops->matinverseadd=SMatInverseAdd;
289 sops->matforwardmultiply=SMatCholeskyForwardMultiply;
290 sops->matseturmat=SMatSetURMat;
291 sops->matfull=SMatFull;
292 sops->matdestroy=SMatDestroy;
293 sops->matgetsize=SMatGetSize;
294 sops->matview=SMatView;
295 sops->matlogdet=SMatLogDet;
296 sops->matname=tmatname;
300static int dcholmatcreate(
int n,
char UPLQ, chfac *sp,
304 DSDPCALLOC1(&S,spmat,&info);DSDPCHKERR(info);
305 S->UPLQ=UPLQ; S->n=n; S->sinv=0; S->dsinv=0; S->spsym=sp;
306 info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info);
312static int dcholmatsinverse(
int n, spmat *S1, spmat *S2){
315 DSDPCALLOC2(&ssinv,
double,n*n,&info);
316 S1->sinv=ssinv; S2->sinv=ssinv; S2->dsinv=1;
321#define __FUNCT__ "DSDPDenseDualMatCreate"
322int DSDPDenseDualMatCreate(
int n,
char UPLQ,
329 info=MchlSetup2(n,&sp); DSDPCHKERR(info);
330 info=dcholmatcreate(n,UPLQ,sp,sops1,smat1);DSDPCHKERR(info);
331 info=MchlSetup2(n,&sp); DSDPCHKERR(info);
332 info=dcholmatcreate(n,UPLQ,sp,sops1,smat2);DSDPCHKERR(info);
333 info=dcholmatsinverse(n,(spmat*)(*smat1),(spmat*)(*smat2));DSDPCHKERR(info);
335 DSDPFunctionReturn(0);
340#define __FUNCT__ "DSDPSparseDualMatCreate"
341int DSDPSparseDualMatCreate(
int n,
int *rnnz,
int *snnz,
342 int trank,
char UPLQ,
int*nnzz,
349 SymbProc(rnnz,snnz,n,&sp); DSDPCHKERR(info);
350 info=dcholmatcreate(n,UPLQ,sp,sops1,smat1);DSDPCHKERR(info);
351 SymbProc(rnnz,snnz,n,&sp); DSDPCHKERR(info);
352 info=dcholmatcreate(n,UPLQ,sp,sops2,smat2);DSDPCHKERR(info);
353 nnz=sp->unnz;*nnzz=nnz;
355 info=dcholmatsinverse(n,(spmat*)(*smat1),(spmat*)(*smat2));DSDPCHKERR(info);
357 DSDPFunctionReturn(0);
int DSDPDualMatOpsInitialize(struct DSDPDualMat_Ops *sops)
Set pointers to null.
Structure of function pointers that each symmetric positive definite matrix type (dense,...
DSDP uses BLAS and LAPACK for many of its operations.
Error handling, printing, and profiling.
Table of function pointers that operate on the S matrix.