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 
11 int DSDPSparsityInSchurMat(DSDP, int, int *, int);
12 int DSDPSetSchurMatOps(DSDP,struct DSDPSchurMat_Ops*,void*);
13 
14 typedef 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 
26 static int dsdpuselapack=1;
27 #undef __FUNCT__
28 #define __FUNCT__ "DSDPUseLAPACKForSchur"
29 int DSDPUseLAPACKForSchur(DSDP dsdp,int flag){
30  DSDPFunctionBegin;
31  dsdpuselapack = flag;
32  DSDPFunctionReturn(0);
33 }
34 
35 #undef __FUNCT__
36 #define __FUNCT__ "Taddline"
37 static 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"
47 static 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"
57 static int Tassemble(void*M){
58  DSDPFunctionBegin;
59  DSDPFunctionReturn(0);
60 }
61 
62 
63 #undef __FUNCT__
64 #define __FUNCT__ "DSDPCheckForSparsity"
65 static 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"
81 static 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"
134 static 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"
154 static 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"
171 static 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"
192 static 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"
202 static 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 
216 static int Tsetup(void*M, int m){
217  DSDPFunctionBegin;
218  DSDPFunctionReturn(0);
219 }
220 
221 static 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 
229 static 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 
237 static 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 
245 static 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 /*
255 static 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 */
263 static const char* tmatname="SPARSE PSD";
264 
265 static 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 
287 static struct DSDPSchurMat_Ops dsdpmatops;
288 
289 int DSDPGetDiagSchurMat(int, struct DSDPSchurMat_Ops **,void **);
290 int DSDPGetLAPACKSUSchurOps(int,struct DSDPSchurMat_Ops**,void**);
291 int DSDPGetLAPACKPUSchurOps(int,struct DSDPSchurMat_Ops**,void**);
292 
293 static int DSDPCreateM(MCholSolverALL*,chfac**,int[],int[],int);
294 static int DSDPCreateSchurMatrix(void*,int);
295 
296 #undef __FUNCT__
297 #define __FUNCT__ "DSDPCreateSchurMatrix"
298 static 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 
360 static struct DSDPSchurMat_Ops dsdpmatops000;
361 
362 #undef __FUNCT__
363 #define __FUNCT__ "DSDPSetDefaultSchurMatrixStructure"
364 int 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
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition: dsdpvec.h:25
Error handling, printing, and profiling.
int DSDPSparsityInSchurMat(DSDP dsdp, int row, int rnnz[], int mm)
Identify nonzero elements in a row of the Schur complement.
Definition: dsdpschurmat.c:649
struct DSDP_C * DSDP
An implementation of the dual-scaling algorithm for semidefinite programming.
Function pointers that a Schur complement matrix (dense, sparse, parallel dense) must provide...
Internal structures for the DSDP solver.
Definition: dsdp.h:65
The API to DSDP for those applications using DSDP as a subroutine library.
Vector operations used by the solver.
int DSDPSchurMatOpsInitialize(struct DSDPSchurMat_Ops *dops)
Initialize function pointers to 0.
Definition: dsdpschurmat.c:44