DSDP
cholmat.c
Go to the documentation of this file.
1#include "numchol.h"
2#include "dsdpschurmat_impl.h"
3#include "dsdpsys.h"
4#include "dsdpvec.h"
5#include "dsdp5.h"
6
11int DSDPSparsityInSchurMat(DSDP, int, int *, int);
12int DSDPSetSchurMatOps(DSDP,struct DSDPSchurMat_Ops*,void*);
13
14typedef struct {
15 chfac *M;
16 int m;
17 int isdense;
18 int *rnnz;
19 int *colnnz;
20 int nnz;
21 DSDPVec D1;
22 DSDP dsdp;
23} MCholSolverALL;
24
25
26static int dsdpuselapack=1;
27#undef __FUNCT__
28#define __FUNCT__ "DSDPUseLAPACKForSchur"
29int DSDPUseLAPACKForSchur(DSDP dsdp,int flag){
30 DSDPFunctionBegin;
31 dsdpuselapack = flag;
32 DSDPFunctionReturn(0);
33}
34
35#undef __FUNCT__
36#define __FUNCT__ "Taddline"
37static int Taddline(void* M, int row, double dd, double v[], int m){
38 int info;
39 MCholSolverALL *AMA = (MCholSolverALL *)M;
40 DSDPFunctionBegin;
41 info=MatAddColumn4(AMA->M,dd,v,row); DSDPCHKERR(info);
42 DSDPFunctionReturn(0);
43}
44
45#undef __FUNCT__
46#define __FUNCT__ "Tadddiagonal"
47static int Tadddiagonal(void* M, int row, double v){
48 int info;
49 MCholSolverALL *AMA = (MCholSolverALL *)M;
50 DSDPFunctionBegin;
51 info=MatAddDiagonalElement(AMA->M,row,v); DSDPCHKERR(info);
52 DSDPFunctionReturn(0);
53}
54
55#undef __FUNCT__
56#define __FUNCT__ "Tassemble"
57static int Tassemble(void*M){
58 DSDPFunctionBegin;
59 DSDPFunctionReturn(0);
60}
61
62
63#undef __FUNCT__
64#define __FUNCT__ "DSDPCheckForSparsity"
65static int DSDPCheckForSparsity( DSDP dsdp, int m, int *rnnz, int tnnz[], int *totalnnz){
66 int info,i,j,tottalnnz=0;
67 DSDPFunctionBegin;
68 memset(rnnz,0,(m+1)*sizeof(int));
69 for (i=0;i<m;i++){
70 info=DSDPSparsityInSchurMat(dsdp,i,tnnz,m); DSDPCHKERR(info);
71 for (j=i+1;j<m;j++){ if (tnnz[j]>0) rnnz[i+1]++;}
72 }
73
74 for (i=0;i<m;i++){tottalnnz+=rnnz[i+1];}
75 *totalnnz=tottalnnz;
76 DSDPFunctionReturn(0);
77}
78
79#undef __FUNCT__
80#define __FUNCT__ "DSDPCreateM"
81static int DSDPCreateM(MCholSolverALL *ABA, chfac **M, int rrnnz[],int tnnz[], int totalnnz){
82
83 int *snnz,*rnnz;
84 int i,j,tt,info;
85 int col,*perm,k;
86 chfac * sfptr;
87 int n=ABA->m,m=ABA->m;
88 DSDP dsdp=ABA->dsdp;
89
90 DSDPFunctionBegin;
91
92 DSDPCALLOC2(&snnz,int,(totalnnz+1),&info); DSDPCHKERR(info);
93 DSDPCALLOC2(&rnnz,int,(m+1),&info); DSDPCHKERR(info);
94 for (i=0;i<=m;i++){ rnnz[i]=rrnnz[i];}
95 tt=0;
96 for (i=0;i<m;i++){
97 info=DSDPSparsityInSchurMat(dsdp,i,tnnz,m); DSDPCHKERR(info);
98 for (j=i+1;j<m;j++){ if (tnnz[j]>0) {snnz[tt]=j; tt++;} }
99 }
100
101 printf("Trying Sparse M: Total nonzeros: %d of %d \n",totalnnz,m*(m-1)/2 );
102 /* Create sparse dual matrix structure */
103 SymbProc(rnnz+1,snnz,n,&sfptr);
104 ABA->isdense=0;
105 ABA->M=sfptr;
106 ABA->nnz=totalnnz;
107 ABA->rnnz=rnnz;
108 ABA->colnnz=snnz;
109 *M=sfptr;
110
111 for (i=0;i<m;i++){
112 rnnz[i+1]+=rnnz[i];
113 }
114
115 perm=sfptr->invp;
116 for (i=m-1;i>=0;i--){
117 for (j=rnnz[i+1]-1;j>=rnnz[i];j--){
118 col=snnz[j];
119 tt=perm[col];
120 if (perm[i] > tt ){
121 for (k=j;k<rnnz[col]-1;k++){ snnz[k]=snnz[k+1];}
122 for (k=i;k<col;k++){ rnnz[k+1]--;}
123 snnz[rnnz[col]]=i;
124 }
125 }
126 }
127
128 DSDPFunctionReturn(0);
129}
130
131
132#undef __FUNCT__
133#define __FUNCT__ "DSDPLinearSolverPrepare"
134static int DSDPLinearSolverPrepare(void* ctx,int*flag){
135
136 cfc_sta Cfact;
137 chfac *sf;
138 MCholSolverALL *AMA = (MCholSolverALL *)ctx;
139
140 DSDPFunctionBegin;
141 *flag=0;
142 sf=AMA->M;
143 /*
144 Cfact=(cfc_sta)ChlFact(sf,sf->iw,sf->rw,FALSE);
145 */
146 Cfact=(cfc_sta)ChlFact(sf,sf->iw,sf->rw,TRUE);
147 if (CfcOk!=Cfact ){ *flag=1; /* printf("Not Good factorization \n"); */ }
148 DSDPFunctionReturn(0);
149}
150
151
152#undef __FUNCT__
153#define __FUNCT__ "DSDPLinearSolve"
154static int DSDPLinearSolve2(void* ctx, double dd[], double xx[], int n){
155
156 int i,info;
157 double *bb;
158 MCholSolverALL *AMA = (MCholSolverALL *)ctx;
159
160 DSDPFunctionBegin;
161 info=DSDPVecGetArray(AMA->D1, &bb);DSDPCHKERR(info);
162 for (i=0;i<n;i++){ bb[i]=dd[i];}
163 ChlSolve(AMA->M, bb, xx);
164 info=DSDPVecRestoreArray(AMA->D1, &bb);
165 DSDPFunctionReturn(0);
166}
167
168
169#undef __FUNCT__
170#define __FUNCT__ "DSDPGramMatRowNonzeros"
171static int DSDPGramMatRowNonzeros(void*M, int row, double cols[], int *ncols, int nrows){
172 MCholSolverALL *AMA = (MCholSolverALL *)M;
173 int i;
174 DSDPFunctionBegin;
175 if (AMA->isdense){
176 *ncols = nrows-row;
177 for (i=row;i<nrows;i++){
178 cols[i]=1.0;
179 }
180 } else {
181 *ncols = AMA->rnnz[row+1] - AMA->rnnz[row]+1;
182 cols[row]=1.0;
183 for (i=AMA->rnnz[row]; i<AMA->rnnz[row+1]; i++){
184 cols[AMA->colnnz[i]]=1.0;
185 }
186 }
187 DSDPFunctionReturn(0);
188}
189
190#undef __FUNCT__
191#define __FUNCT__ "Tzero"
192static int Tzero(void*M){
193 int info;
194 MCholSolverALL *AMA = (MCholSolverALL *)M;
195 DSDPFunctionBegin;
196 info=MatZeroEntries4(AMA->M); DSDPCHKERR(info);
197 DSDPFunctionReturn(0);
198}
199
200#undef __FUNCT__
201#define __FUNCT__ "Tdestroy"
202static int Tdestroy(void*M){
203 MCholSolverALL *AMA = (MCholSolverALL *)M;
204 int info;
205 DSDPFunctionBegin;
206 CfcFree(&AMA->M);
207 info = DSDPVecDestroy(&AMA->D1); DSDPCHKERR(info);
208 if (AMA->isdense==0 && AMA->rnnz ){
209 DSDPFREE(&AMA->rnnz,&info);DSDPCHKERR(info);
210 DSDPFREE(&AMA->colnnz,&info);DSDPCHKERR(info);
211 }
212 DSDPFREE(&AMA,&info);DSDPCHKERR(info);
213 DSDPFunctionReturn(0);
214}
215
216static int Tsetup(void*M, int m){
217 DSDPFunctionBegin;
218 DSDPFunctionReturn(0);
219}
220
221static int TTTMatMult(void*M,double x[],double y[],int n){
222 MCholSolverALL *AMA = (MCholSolverALL *)M;
223 int info;
224 DSDPFunctionBegin;
225 info=MatMult4(AMA->M,x,y,n); DSDPCHKERR(info);
226 DSDPFunctionReturn(0);
227}
228
229static int TTTMatShiftDiagonal(void *M, double d){
230 MCholSolverALL *AMA = (MCholSolverALL *)M;
231 int info;
232 DSDPFunctionBegin;
233 info=Mat4DiagonalShift(AMA->M,d); DSDPCHKERR(info);
234 DSDPFunctionReturn(0);
235}
236
237static int TTTMatAddDiagonal(void *M, double diag[], int m){
238 MCholSolverALL *AMA = (MCholSolverALL *)M;
239 int info;
240 DSDPFunctionBegin;
241 info=Mat4AddDiagonal(AMA->M,diag,m); DSDPCHKERR(info);
242 DSDPFunctionReturn(0);
243}
244
245static int TTTMatView(void *M){
246 MCholSolverALL *AMA = (MCholSolverALL *)M;
247 int info;
248 DSDPFunctionBegin;
249 info=Mat4View(AMA->M); DSDPCHKERR(info);
250 DSDPFunctionReturn(0);
251}
252
253
254/*
255static int TTTMatGetDiagonal(void *M, double d[], int n){
256 chfac*A = (chfac*)M;
257 int info;
258 DSDPFunctionBegin;
259 info=Mat4GetDiagonal(A,d,n); DSDPCHKERR(info);
260 DSDPFunctionReturn(0);
261}
262*/
263static const char* tmatname="SPARSE PSD";
264
265static int TMatOpsInit(struct DSDPSchurMat_Ops *mops){
266 int info;
267 DSDPFunctionBegin;
268 info=DSDPSchurMatOpsInitialize(mops); DSDPCHKERR(info);
269 mops->matrownonzeros=DSDPGramMatRowNonzeros;
270 mops->mataddrow=Taddline;
271 mops->matadddiagonal=TTTMatAddDiagonal;
272 mops->mataddelement=Tadddiagonal;
273 mops->matshiftdiagonal=TTTMatShiftDiagonal;
274 mops->matassemble=Tassemble;
275 mops->matscaledmultiply=TTTMatMult;
276 mops->matfactor=DSDPLinearSolverPrepare;
277 mops->matsolve=DSDPLinearSolve2;
278 mops->matdestroy=Tdestroy;
279 mops->matzero=Tzero;
280 mops->matsetup=Tsetup;
281 mops->matview=TTTMatView;
282 mops->id=5;
283 mops->matname=tmatname;
284 DSDPFunctionReturn(0);
285}
286
287static struct DSDPSchurMat_Ops dsdpmatops;
288
289int DSDPGetDiagSchurMat(int, struct DSDPSchurMat_Ops **,void **);
290int DSDPGetLAPACKSUSchurOps(int,struct DSDPSchurMat_Ops**,void**);
291int DSDPGetLAPACKPUSchurOps(int,struct DSDPSchurMat_Ops**,void**);
292
293static int DSDPCreateM(MCholSolverALL*,chfac**,int[],int[],int);
294static int DSDPCreateSchurMatrix(void*,int);
295
296#undef __FUNCT__
297#define __FUNCT__ "DSDPCreateSchurMatrix"
298static int DSDPCreateSchurMatrix(void *ctx, int m){
299
300 int info;
301 int *rnnz,*tnnz,totalnnz;
302 int gotit=0;
303 DSDP dsdp=(DSDP)ctx;
304 chfac *sf;
305 MCholSolverALL *AMA;
306 void *tdata;
307 struct DSDPSchurMat_Ops *tmatops;
308
309 DSDPFunctionBegin;
310 if (m<=1){
311 info=DSDPGetDiagSchurMat(m,&tmatops,&tdata); DSDPCHKERR(info);
312 info=DSDPSetSchurMatOps(dsdp,tmatops,tdata); DSDPCHKERR(info);
313 return 0;
314 }
315
316 DSDPCALLOC2(&rnnz,int,(m+1),&info); DSDPCHKERR(info);
317 DSDPCALLOC2(&tnnz,int,(m+1),&info); DSDPCHKERR(info);
318
319 info=DSDPCheckForSparsity(dsdp,m,rnnz,tnnz,&totalnnz); DSDPCHKERR(info);
320
321 if (totalnnz*2+m > m*m*0.1 && dsdpuselapack) {
322 gotit=1;
323 info=DSDPGetLAPACKSUSchurOps(m,&tmatops,&tdata);
324 /* DSDPCHKERR(info); */ if (info) {gotit=0;printf("Try packed format\n"); }
325 DSDPLogInfo(0,8,"Creating Dense Full LAPACK Schur Matrix\n");
326 info=DSDPSetSchurMatOps(dsdp,tmatops,tdata); DSDPCHKERR(info);
327 }
328
329 if ( 0 && totalnnz*2+m > m*m*0.1 && dsdpuselapack) {
330
331 info=DSDPGetLAPACKPUSchurOps(m,&tmatops,&tdata); DSDPCHKERR(info);
332 DSDPLogInfo(0,8,"Creating Dense Packed LAPACK Schur Matrix\n");
333 info=DSDPSetSchurMatOps(dsdp,tmatops,tdata); DSDPCHKERR(info);
334 gotit=1;
335
336 }
337 if (gotit==0){
338
339 DSDPCALLOC1(&AMA,MCholSolverALL,&info);DSDPCHKERR(info);
340 AMA->dsdp=dsdp; AMA->m=m;
341 info=DSDPVecCreateSeq(m,&AMA->D1); DSDPCHKERR(info);
342 if (totalnnz*2+m > m*m * 0.11 ){
343 info=MchlSetup2(m,&sf); DSDPCHKERR(info);
344 AMA->M=sf; AMA->isdense=1;
345 AMA->rnnz=0; AMA->colnnz=0;
346 DSDPLogInfo(0,8,"Creating Dense Full non LAPACK Schur Matrix\n");
347 } else {
348 info=DSDPCreateM(AMA,&sf,rnnz,tnnz,totalnnz); DSDPCHKERR(info);
349 DSDPLogInfo(0,8,"Creating Sparse Schur Matrix\n");
350 }
351 AMA->M=sf;
352 info=TMatOpsInit(&dsdpmatops); DSDPCHKERR(info);
353 info=DSDPSetSchurMatOps(dsdp,&dsdpmatops,(void*)AMA); DSDPCHKERR(info);
354 }
355 DSDPFREE(&tnnz,&info);DSDPCHKERR(info);
356 DSDPFREE(&rnnz,&info);DSDPCHKERR(info);
357 DSDPFunctionReturn(0);
358}
359
360static struct DSDPSchurMat_Ops dsdpmatops000;
361
362#undef __FUNCT__
363#define __FUNCT__ "DSDPSetDefaultSchurMatrixStructure"
364int DSDPSetDefaultSchurMatrixStructure(DSDP dsdp){
365 int info;
366 DSDPFunctionBegin;
367 info=DSDPSchurMatOpsInitialize(&dsdpmatops000); DSDPCHKERR(info);
368 dsdpmatops000.matsetup=DSDPCreateSchurMatrix;
369 info=DSDPSetSchurMatOps(dsdp,&dsdpmatops000,(void*)dsdp);DSDPCHKERR(info);
370 DSDPFunctionReturn(0);
371}
372
int DSDPSetSchurMatOps(DSDP, struct DSDPSchurMat_Ops *, void *)
Set the Schur complement matrix.
Definition dsdpcops.c:602
The API to DSDP for those applications using DSDP as a subroutine library.
struct DSDP_C * DSDP
An implementation of the dual-scaling algorithm for semidefinite programming.
int DSDPSparsityInSchurMat(DSDP dsdp, int row, int rnnz[], int mm)
Identify nonzero elements in a row of the Schur complement.
int DSDPSchurMatOpsInitialize(struct DSDPSchurMat_Ops *dops)
Initialize function pointers to 0.
Function pointers that a Schur complement matrix (dense, sparse, parallel dense) must provide.
Error handling, printing, and profiling.
Vector operations used by the solver.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition dsdpvec.h:25
Internal structures for the DSDP solver.
Definition dsdp.h:65