00001
00002 #include "dsdpschurmat_impl.h"
00003 #include "dsdpdualmat_impl.h"
00004 #include "dsdpdsmat_impl.h"
00005 #include "dsdpsys.h"
00012 typedef struct {
00013 int n;
00014 double *val;
00015 int owndata;
00016 } diagmat;
00017
00018 static int DiagMatCreate(int,diagmat**);
00019 static int DiagMatMult(void*,double[],double[],int);
00020 static int DiagMatGetSize(void*, int *);
00021 static int DiagMatAddRow2(void*, int, double, double[], int);
00022 static int DiagMatDestroy(void*);
00023 static int DiagMatView(void*);
00024 static int DiagMatLogDeterminant(void*, double *);
00025
00026
00027
00028 static int DiagMatCreate(int n,diagmat**M){
00029 int info;
00030 diagmat *M7;
00031
00032 DSDPCALLOC1(&M7,diagmat,&info);DSDPCHKERR(info);
00033 DSDPCALLOC2(&M7->val,double,n,&info);DSDPCHKERR(info);
00034 if (n>0 && M7->val==NULL) return 1;
00035 M7->owndata=1; M7->n=n;
00036 *M=M7;
00037 return 0;
00038 }
00039
00040 static int DiagMatDestroy(void* AA){
00041 int info;
00042 diagmat* A=(diagmat*) AA;
00043 if (A->owndata && A->val){ DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
00044 DSDPFREE(&A,&info);DSDPCHKERR(info);
00045 return 0;
00046 }
00047
00048
00049 static int DiagMatMult(void* AA, double x[], double y[], int n){
00050
00051 diagmat* A=(diagmat*) AA;
00052 double *vv=A->val;
00053 int i;
00054
00055 if (A->n != n) return 1;
00056 if (x==0 && n>0) return 3;
00057 if (y==0 && n>0) return 3;
00058
00059 for (i=0; i<n; i++){
00060 y[i]=x[i]*vv[i];
00061 }
00062 return 0;
00063 }
00064
00065
00066 static int DiagMatGetSize(void *AA, int *n){
00067 diagmat* A=(diagmat*) AA;
00068 *n=A->n;
00069 return 0;
00070 }
00071
00072 static int DiagMatView(void* AA){
00073 diagmat* A=(diagmat*) AA;
00074 int i;
00075 for (i=0;i<A->n; i++){
00076 printf(" Row: %d, Column: %d, Value: %8.4e \n",i,i,A->val[i]);
00077 }
00078 return 0;
00079 }
00080
00081 static int DiagMatAddRow2(void* AA, int nrow, double dd, double row[], int n){
00082 diagmat* A=(diagmat*) AA;
00083 A->val[nrow] += dd*row[nrow];
00084 return 0;
00085 }
00086
00087
00088 static int DiagMatAddElement(void*A, int nrow, double dd){
00089 diagmat* AA = (diagmat*)A;
00090 AA->val[nrow] += dd;
00091 return 0;
00092 }
00093
00094 static int DiagMatZeros(void*A){
00095 diagmat* AA = (diagmat*)A;
00096 int n=AA->n;
00097 memset(AA->val,0,n*sizeof(double));
00098 return 0;
00099 }
00100
00101 static int DiagMatSolve(void* A, double b[], double x[],int n){
00102 diagmat* AA = (diagmat*)A;
00103 double *v=AA->val;
00104 int i;
00105 for (i=0;i<n;i++){
00106 x[i]=b[i]/v[i];
00107 }
00108 return 0;
00109 }
00110
00111 static int DiagMatSolve2(void* A, int indx[], int nindx, double b[], double x[],int n){
00112 diagmat* AA = (diagmat*)A;
00113 double *v=AA->val;
00114 int i,j;
00115 memset((void*)x,0,n*sizeof(double));
00116 for (j=0;j<nindx;j++){
00117 i=indx[j];
00118 x[i]=b[i]/v[i];
00119 }
00120 return 0;
00121 }
00122
00123 static int DiagMatCholeskySolveBackward(void* A, double b[], double x[],int n){
00124 int i;
00125 for (i=0;i<n;i++){
00126 x[i]=b[i];
00127 }
00128 return 0;
00129 }
00130
00131 static int DiagMatInvert(void *A){
00132 return 0;
00133 }
00134
00135 static int DiagMatCholeskyFactor(void*A,int *flag){
00136 diagmat* AA = (diagmat*)A;
00137 double *v=AA->val;
00138 int i,n=AA->n;
00139 *flag=0;
00140 for (i=0;i<n;i++){
00141 if (v[i]<=0){ *flag=i+1; break;}
00142 }
00143 return 0;
00144 }
00145
00146 static int DiagMatLogDeterminant(void*A, double *dd){
00147 diagmat* AA = (diagmat*)A;
00148 double d=0,*val=AA->val;
00149 int i,n=AA->n;
00150 for (i=0;i<n;i++){
00151 if (val[i]<=0) return 1;
00152 d+=log(val[i]);
00153 }
00154 *dd=d;
00155 return 0;
00156 }
00157
00158 static int DiagMatTakeUREntriesP(void*A, double dd[], int nn, int n){
00159 diagmat* AA = (diagmat*)A;
00160 int i,ii;
00161 double *val=AA->val;
00162 for (i=0;i<n;i++){
00163 ii=(i+1)*(i+2)/2-1;
00164 val[i]=dd[ii];
00165 }
00166 return 0;
00167 }
00168 static int DiagMatTakeUREntriesU(void*A, double dd[], int nn, int n){
00169 diagmat* AA = (diagmat*)A;
00170 int i;
00171 double *val=AA->val;
00172 for (i=0;i<n;i++){
00173 val[i]=dd[i*n+i];
00174 }
00175 return 0;
00176 }
00177
00178 static int DiagMatInverseAddP(void*A, double alpha, double dd[], int nn, int n){
00179 diagmat* AA = (diagmat*)A;
00180 int i,ii;
00181 double *val=AA->val;
00182 for (i=0;i<n;i++){
00183 ii=(i+1)*(i+2)/2-1;
00184 dd[ii]+=alpha/val[i];
00185 }
00186 return 0;
00187 }
00188 static int DiagMatInverseAddU(void*A, double alpha, double dd[], int nn, int n){
00189 diagmat* AA = (diagmat*)A;
00190 int i;
00191 double *val=AA->val;
00192 for (i=0;i<n;i++){
00193 dd[i*n+i]+=alpha/val[i];
00194 }
00195 return 0;
00196 }
00197
00198 static int DiagMatFull(void*A, int* dfull){
00199 *dfull=1;
00200 return 0;
00201 }
00202
00203 static struct DSDPDualMat_Ops sdmatopsp;
00204 static struct DSDPDualMat_Ops sdmatopsu;
00205 static const char* diagmatname="DIAGONAL";
00206
00207 static int DiagDualOpsInitializeP(struct DSDPDualMat_Ops* sops){
00208 int info;
00209 if (sops==NULL) return 0;
00210 info=DSDPDualMatOpsInitialize(sops); DSDPCHKERR(info);
00211 sops->matcholesky=DiagMatCholeskyFactor;
00212 sops->matsolveforward=DiagMatSolve;
00213 sops->matsolvebackward=DiagMatCholeskySolveBackward;
00214 sops->matinvert=DiagMatInvert;
00215 sops->matinverseadd=DiagMatInverseAddP;
00216 sops->matinversemultiply=DiagMatSolve2;
00217 sops->matseturmat=DiagMatTakeUREntriesP;
00218 sops->matfull=DiagMatFull;
00219 sops->matdestroy=DiagMatDestroy;
00220 sops->matgetsize=DiagMatGetSize;
00221 sops->matview=DiagMatView;
00222 sops->matlogdet=DiagMatLogDeterminant;
00223 sops->id=9;
00224 sops->matname=diagmatname;
00225 return 0;
00226 }
00227 static int DiagDualOpsInitializeU(struct DSDPDualMat_Ops* sops){
00228 int info;
00229 if (sops==NULL) return 0;
00230 info=DSDPDualMatOpsInitialize(sops); DSDPCHKERR(info);
00231 sops->matcholesky=DiagMatCholeskyFactor;
00232 sops->matsolveforward=DiagMatSolve;
00233 sops->matsolvebackward=DiagMatCholeskySolveBackward;
00234 sops->matinvert=DiagMatInvert;
00235 sops->matinversemultiply=DiagMatSolve2;
00236 sops->matseturmat=DiagMatTakeUREntriesU;
00237 sops->matfull=DiagMatFull;
00238 sops->matinverseadd=DiagMatInverseAddU;
00239 sops->matdestroy=DiagMatDestroy;
00240 sops->matgetsize=DiagMatGetSize;
00241 sops->matview=DiagMatView;
00242 sops->matlogdet=DiagMatLogDeterminant;
00243 sops->id=9;
00244 sops->matname=diagmatname;
00245 return 0;
00246 }
00247
00248 #undef __FUNCT__
00249 #define __FUNCT__ "DSDPDiagDualMatCreateP"
00250 int DSDPDiagDualMatCreateP(int n,
00251 struct DSDPDualMat_Ops **sops1, void**smat1,
00252 struct DSDPDualMat_Ops **sops2, void**smat2){
00253 diagmat *AA;
00254 int info;
00255 DSDPFunctionBegin;
00256
00257 info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00258 info=DiagDualOpsInitializeP(&sdmatopsp); DSDPCHKERR(info);
00259 *sops1=&sdmatopsp;
00260 *smat1=(void*)AA;
00261
00262 info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00263 info=DiagDualOpsInitializeP(&sdmatopsp); DSDPCHKERR(info);
00264 *sops2=&sdmatopsp;
00265 *smat2=(void*)AA;
00266 DSDPFunctionReturn(0);
00267 }
00268
00269 #undef __FUNCT__
00270 #define __FUNCT__ "DSDPDiagDualMatCreateU"
00271 int DSDPDiagDualMatCreateU(int n,
00272 struct DSDPDualMat_Ops **sops1, void**smat1,
00273 struct DSDPDualMat_Ops **sops2, void**smat2){
00274 diagmat *AA;
00275 int info;
00276 DSDPFunctionBegin;
00277 info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00278 info=DiagDualOpsInitializeU(&sdmatopsu); DSDPCHKERR(info);
00279 *sops1=&sdmatopsu;
00280 *smat1=(void*)AA;
00281 info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00282 info=DiagDualOpsInitializeU(&sdmatopsu); DSDPCHKERR(info);
00283 *sops2=&sdmatopsu;
00284 *smat2=(void*)AA;
00285 DSDPFunctionReturn(0);
00286 }
00287
00288
00289
00290 static int DiagMatVecVec(void*A,double x[], int n, double *vv){
00291 diagmat* AA = (diagmat*)A;
00292 double *v=AA->val,vAv=0;
00293 int i;
00294 for (i=0;i<n;i++){
00295 vAv+=x[i]*x[i]*v[i];
00296 }
00297 *vv=vAv;
00298 return 0;
00299 }
00300
00301 static int DDiagDSMatOpsInitP(struct DSDPDSMat_Ops *ddiagops){
00302 int info;
00303 if (ddiagops==NULL) return 0;
00304 info=DSDPDSMatOpsInitialize(ddiagops);DSDPCHKERR(info);
00305 ddiagops->matseturmat=DiagMatTakeUREntriesP;
00306 ddiagops->matview=DiagMatView;
00307 ddiagops->matgetsize=DiagMatGetSize;
00308 ddiagops->matmult=DiagMatMult;
00309 ddiagops->matvecvec=DiagMatVecVec;
00310 ddiagops->matzeroentries=DiagMatZeros;
00311 ddiagops->matdestroy=DiagMatDestroy;
00312 ddiagops->id=9;
00313 ddiagops->matname=diagmatname;
00314 DSDPFunctionReturn(0);
00315 }
00316 static int DDiagDSMatOpsInitU(struct DSDPDSMat_Ops *ddiagops){
00317 int info;
00318 if (ddiagops==NULL) return 0;
00319 info=DSDPDSMatOpsInitialize(ddiagops);DSDPCHKERR(info);
00320 ddiagops->matseturmat=DiagMatTakeUREntriesU;
00321 ddiagops->matview=DiagMatView;
00322 ddiagops->matgetsize=DiagMatGetSize;
00323 ddiagops->matmult=DiagMatMult;
00324 ddiagops->matvecvec=DiagMatVecVec;
00325 ddiagops->matzeroentries=DiagMatZeros;
00326 ddiagops->matdestroy=DiagMatDestroy;
00327 ddiagops->id=9;
00328 ddiagops->matname=diagmatname;
00329 DSDPFunctionReturn(0);
00330 }
00331
00332 static struct DSDPDSMat_Ops dsdiagmatopsp;
00333 static struct DSDPDSMat_Ops dsdiagmatopsu;
00334
00335 #undef __FUNCT__
00336 #define __FUNCT__ "DSDPDiagDSMatP"
00337 int DSDPCreateDiagDSMatP(int n,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
00338
00339 int info=0;
00340 diagmat *AA;
00341
00342 DSDPFunctionBegin;
00343 info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00344 info=DDiagDSMatOpsInitP(&dsdiagmatopsp); DSDPCHKERR(info);
00345 *dsmatops=&dsdiagmatopsp;
00346 *dsmat=(void*)AA;
00347 DSDPFunctionReturn(0);
00348 }
00349 #undef __FUNCT__
00350 #define __FUNCT__ "DSDPDiagDSMatU"
00351 int DSDPCreateDiagDSMatU(int n,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
00352
00353 int info=0;
00354 diagmat *AA;
00355
00356 DSDPFunctionBegin;
00357 info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00358 info=DDiagDSMatOpsInitU(&dsdiagmatopsu); DSDPCHKERR(info);
00359 *dsmatops=&dsdiagmatopsu;
00360 *dsmat=(void*)AA;
00361 DSDPFunctionReturn(0);
00362 }
00363
00364 static int DiagRowNonzeros(void*M, int row, double cols[], int *ncols,int nrows){
00365 DSDPFunctionBegin;
00366 *ncols = 1;
00367 cols[row]=1.0;
00368 DSDPFunctionReturn(0);
00369 }
00370
00371 static int DiagAddDiag(void*M, double diag[], int m){
00372 diagmat* AA = (diagmat*)M;
00373 double *v=AA->val;
00374 int i;
00375 DSDPFunctionBegin;
00376 for (i=0;i<m;i++){
00377 v[i]+=diag[i];
00378 }
00379 DSDPFunctionReturn(0);
00380 }
00381
00382 static int DiagMultiply(void*M, double xin[], double xout[], int m){
00383 diagmat* AA = (diagmat*)M;
00384 double *v=AA->val;
00385 int i;
00386 DSDPFunctionBegin;
00387 for (i=0;i<m;i++){
00388 xout[i]+=v[i]*xin[i];
00389 }
00390 DSDPFunctionReturn(0);
00391 }
00392
00393 static int DiagShiftDiag(void*M, double dd){
00394 diagmat* AA = (diagmat*)M;
00395 double *v=AA->val;
00396 int i,m=AA->n;
00397 DSDPFunctionBegin;
00398 for (i=0;i<m;i++){
00399 v[i]+=dd;
00400 }
00401 DSDPFunctionReturn(0);
00402 }
00403
00404 static int DiagAddElement(void*M, int ii, double dd){
00405 diagmat* AA = (diagmat*)M;
00406 DSDPFunctionBegin;
00407 AA->val[ii]+=dd;
00408 DSDPFunctionReturn(0);
00409 }
00410
00411 static int DiagMatOnProcessor(void*A,int row,int*flag){
00412 *flag=1;
00413 return 0;
00414 }
00415
00416 static int DiagAssemble(void*M){
00417 return 0;
00418 }
00419
00420 static struct DSDPSchurMat_Ops dsdpdiagschurops;
00421
00422 #undef __FUNCT__
00423 #define __FUNCT__ "DSDPDiagSchurOps"
00424 static int DiagSchurOps(struct DSDPSchurMat_Ops *sops){
00425 int info;
00426 DSDPFunctionBegin;
00427 if (!sops) return 0;
00428 info=DSDPSchurMatOpsInitialize(sops); DSDPCHKERR(info);
00429 sops->matzero=DiagMatZeros;
00430 sops->mataddrow=DiagMatAddRow2;
00431 sops->mataddelement=DiagMatAddElement;
00432 sops->matdestroy=DiagMatDestroy;
00433 sops->matfactor=DiagMatCholeskyFactor;
00434 sops->matsolve=DiagMatSolve;
00435 sops->matadddiagonal=DiagAddDiag;
00436 sops->matshiftdiagonal=DiagShiftDiag;
00437 sops->mataddelement=DiagAddElement;
00438 sops->matscaledmultiply=DiagMultiply;
00439 sops->matassemble=DiagAssemble;
00440 sops->pmatonprocessor=DiagMatOnProcessor;
00441 sops->matrownonzeros=DiagRowNonzeros;
00442 sops->id=9;
00443 sops->matname=diagmatname;
00444 DSDPFunctionReturn(0);
00445 }
00446
00447 #undef __FUNCT__
00448 #define __FUNCT__ "DSDPGetDiagSchurMat"
00449 int DSDPGetDiagSchurMat(int n, struct DSDPSchurMat_Ops **sops, void **data){
00450 int info=0;
00451 diagmat *AA;
00452 DSDPFunctionBegin;
00453 info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00454 info=DiagSchurOps(&dsdpdiagschurops); DSDPCHKERR(info);
00455 if (sops){*sops=&dsdpdiagschurops;}
00456 if (data){*data=(void*)AA;}
00457 DSDPFunctionReturn(0);
00458 }