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.
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.