20 static int SMatDestroy(
void*S){
25 DSDPFREE(&SS->sinv,&info);
31 static int SMatGetSize(
void *S,
int *n){
37 static int SMatView(
void* S){
40 info=Mat4View(SS->spsym); DSDPCHKERR(info);
44 static int SMatLogDet(
void* S,
double *dd){
47 info=Mat4LogDet(SS->spsym,dd); DSDPCHKERR(info);
53 static 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);
92 static 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);
129 static 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);
140 static 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);
159 static int SMatCholeskySolveBackward(
void *S,
double b[],
double x[],
int n){
161 ChlSolveBackward(SS->spsym, b, x);
165 static int SMatCholeskySolveForward(
void *S,
double b[],
double x[],
int n){
167 ChlSolveForward(SS->spsym, b, x);
171 static int SMatFull(
void *S,
int *full){
176 static int SMatCholeskyFactor(
void *S,
int *flag){
183 Cfact=(cfc_sta)ChlFact(SS->spsym,iw,rw,TRUE);
192 static 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);
219 static 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);
242 static 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);
253 static int SMatCholeskyForwardMultiply(
void *S,
double b[],
double x[],
int n){
255 GetUhat(SS->spsym,b,x);
259 static 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));
278 static 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;
300 static 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);
312 static 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" 322 int 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" 341 int 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.
Error handling, printing, and profiling.
Table of function pointers that operate on the S matrix.
Structure of function pointers that each symmetric positive definite matrix type (dense, sparse) must implement.
DSDP uses BLAS and LAPACK for many of its operations.