DSDP
|
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 /* static int DiagMatScale(double *, diagmat *); */ 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 }