DSDP
|
00001 #include "dsdpsys.h" 00002 #include "dsdpdsmat_impl.h" 00003 00008 typedef struct { 00009 int n; 00010 double *an; 00011 int *col; 00012 int *nnz; 00013 } spdsmat; 00014 00015 static int SpSymMatSetURValuesP(void*DS, double v[], int nn, int n){ 00016 spdsmat*ds=(spdsmat*)DS; 00017 int i,j,k1,k2,*nnz=ds->nnz,*col=ds->col; 00018 double *an=ds->an; 00019 for (i=0;i<n;i++,nnz++){ 00020 k1=*nnz; k2=*(nnz+1); 00021 for (j=k1;j<k2;j++,an++,col++){ 00022 if ((*col)==i){ *an = v[*col]/2;} 00023 else { *an = v[*col]; } 00024 } 00025 v+=i+1; 00026 } 00027 return 0; 00028 } 00029 00030 static int SpSymMatSetURValuesU(void*DS, double v[], int nn, int n){ 00031 spdsmat*ds=(spdsmat*)DS; 00032 int i,j,k1,k2,*nnz=ds->nnz,*col=ds->col; 00033 double *an=ds->an; 00034 for (i=0;i<n;i++,nnz++){ 00035 k1=*nnz; k2=*(nnz+1); 00036 for (j=k1;j<k2;j++,an++,col++){ 00037 if ((*col)==i){ *an = v[*col]/2;} 00038 else { *an = v[*col]; } 00039 } 00040 v+=n; 00041 } 00042 return 0; 00043 } 00044 00045 static int SpSymMatView(void *DS){ 00046 spdsmat*ds=(spdsmat*)DS; 00047 int i,j,k1,k2,n=ds->n,*nnz=ds->nnz,*col=ds->col; 00048 double *an=ds->an; 00049 for (i=0;i<n;i++){ 00050 k1=nnz[i]; k2=nnz[i+1]; 00051 printf("Row %d: ",i); 00052 for (j=k1;j<k2;j++){ 00053 if (col[j]==i){ printf("%d: %4.4f",col[j],2*an[j]); } 00054 else { printf("%d: %4.4f",col[j],an[j]);} 00055 } 00056 printf("\n"); 00057 } 00058 return 0; 00059 } 00060 /* 00061 static int SpSymMatShiftDiagonal(void *DS, double dd){ 00062 spdsmat*ds=(spdsmat*)DS; 00063 int i,n=ds->n,*nnz=ds->nnz; 00064 double *an=ds->an; 00065 for (i=0;i<n;i++){ 00066 an[nnz[i+1]-1] += dd/2; 00067 } 00068 return 0; 00069 } 00070 */ 00071 static int SpSymMatDestroy(void *DS){ 00072 spdsmat*ds=(spdsmat*)DS; 00073 int info; 00074 DSDPFREE(&ds->nnz,&info);if (info) return 1; 00075 DSDPFREE(&ds->col,&info);if (info) return 1; 00076 DSDPFREE(&ds->an,&info);if (info) return 1; 00077 DSDPFREE(&ds,&info);if (info) return 1; 00078 return 0; 00079 } 00080 00081 static int SpSymMatGetSize(void *DS, int*n){ 00082 spdsmat*ds=(spdsmat*)DS; 00083 *n=ds->n; 00084 return 0; 00085 } 00086 00087 static int SpSymMatZero(void*DS){ 00088 spdsmat*ds=(spdsmat*)DS; 00089 int nn=ds->nnz[ds->n]; 00090 double *an=ds->an; 00091 memset((void*)an,0,nn*sizeof(double)); 00092 return 0; 00093 } 00094 00095 static int SpSymMatMult(void*DS, double x[], double y[], int n){ 00096 spdsmat*ds=(spdsmat*)DS; 00097 int i,j,k1,k2,*nnz=ds->nnz,*col=ds->col; 00098 double *an=ds->an; 00099 memset((void*)y,0,n*sizeof(double)); 00100 for (i=0;i<n;i++,nnz++){ 00101 k1=*nnz; k2=*(nnz+1); 00102 for (j=k1;j<k2;j++,col++,an++){ 00103 y[*col] += x[i] * (*an); 00104 y[i] += x[*col] * (*an); 00105 } 00106 } 00107 return 0; 00108 } 00109 00110 static int SpSymMatVecVec(void*DS, double x[], int n, double *vAv){ 00111 spdsmat*ds=(spdsmat*)DS; 00112 int i,j,k1,k2,*nnz=ds->nnz,*col=ds->col; 00113 double vv,*an=ds->an; 00114 *vAv=0; 00115 for (i=0;i<n;i++,nnz++){ 00116 k1=*nnz; k2=*(nnz+1); 00117 vv=0; 00118 for (j=k1;j<k2;j++,col++,an++){ 00119 vv+=x[*col]*(*an); 00120 } 00121 *vAv+=vv*x[i]*2; 00122 } 00123 return 0; 00124 } 00125 /* 00126 static int SpSymMatAddRow(void *DS, int row, double dd, double v[], int n){ 00127 spdsmat*ds=(spdsmat*)DS; 00128 int j,k1,k2,*nnz=ds->nnz,*col=ds->col; 00129 double *an=ds->an; 00130 k1=nnz[row]; k2=nnz[row+1]; 00131 for (j=k1;j<k2;j++){ 00132 if (row==col[j]){ an[j] += dd*v[col[j]]/2; } 00133 else { an[j] += dd*v[col[j]]; } 00134 } 00135 return 0; 00136 } 00137 */ 00138 static const char* dsmatname="SPARSE, SYMMETRIC MATRIX"; 00139 static int DSDPDSSparseInitializeOpsP(struct DSDPDSMat_Ops* dsops){ 00140 int info; 00141 if (!dsops) return 0; 00142 info=DSDPDSMatOpsInitialize(dsops); DSDPCHKERR(info); 00143 dsops->matseturmat=SpSymMatSetURValuesP; 00144 dsops->matview=SpSymMatView; 00145 dsops->matdestroy=SpSymMatDestroy; 00146 dsops->matgetsize=SpSymMatGetSize; 00147 dsops->matzeroentries=SpSymMatZero; 00148 dsops->matmult=SpSymMatMult; 00149 dsops->matvecvec=SpSymMatVecVec; 00150 dsops->id=6; 00151 dsops->matname=dsmatname; 00152 return 0; 00153 } 00154 static int DSDPDSSparseInitializeOpsU(struct DSDPDSMat_Ops* dsops){ 00155 int info; 00156 if (!dsops) return 0; 00157 info=DSDPDSMatOpsInitialize(dsops); DSDPCHKERR(info); 00158 dsops->matseturmat=SpSymMatSetURValuesU; 00159 dsops->matview=SpSymMatView; 00160 dsops->matdestroy=SpSymMatDestroy; 00161 dsops->matgetsize=SpSymMatGetSize; 00162 dsops->matzeroentries=SpSymMatZero; 00163 dsops->matmult=SpSymMatMult; 00164 dsops->matvecvec=SpSymMatVecVec; 00165 dsops->id=6; 00166 dsops->matname=dsmatname; 00167 return 0; 00168 } 00169 00170 static struct DSDPDSMat_Ops tdsdsopsp; 00171 static struct DSDPDSMat_Ops tdsdsopsu; 00172 #undef __FUNCT__ 00173 #define __FUNCT__ "DSDPCreateSparseDSMat" 00174 int DSDPSparseMatCreatePattern2P(int n, int rnnz[], int cols[], int tnnz,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){ 00175 int i,info; 00176 spdsmat*ds; 00177 DSDPFunctionBegin; 00178 DSDPCALLOC1(&ds,spdsmat,&info);DSDPCHKERR(info); 00179 DSDPCALLOC2(&ds->nnz,int,(n+1),&info);DSDPCHKERR(info); 00180 ds->nnz[0]=0; 00181 for (i=0;i<n;i++) ds->nnz[i+1]=ds->nnz[i]+rnnz[i]; 00182 DSDPCALLOC2(&ds->col,int,tnnz,&info);DSDPCHKERR(info); 00183 DSDPCALLOC2(&ds->an,double,tnnz,&info);DSDPCHKERR(info); 00184 for (i=0;i<tnnz;i++) ds->col[i]=cols[i]; 00185 info=DSDPDSSparseInitializeOpsP(&tdsdsopsp); DSDPCHKERR(info); 00186 *dsmatops=&tdsdsopsp; 00187 *dsmat=(void*)ds; 00188 DSDPFunctionReturn(0); 00189 } 00190 00191 #undef __FUNCT__ 00192 #define __FUNCT__ "DSDPCreateSparseDSMatU" 00193 int DSDPSparseMatCreatePattern2U(int n, int rnnz[], int cols[], int tnnz,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){ 00194 int i,info; 00195 spdsmat*ds; 00196 DSDPFunctionBegin; 00197 DSDPCALLOC1(&ds,spdsmat,&info);DSDPCHKERR(info); 00198 DSDPCALLOC2(&ds->nnz,int,(n+1),&info);DSDPCHKERR(info); 00199 ds->nnz[0]=0; 00200 for (i=0;i<n;i++) ds->nnz[i+1]=ds->nnz[i]+rnnz[i]; 00201 DSDPCALLOC2(&ds->col,int,tnnz,&info);DSDPCHKERR(info); 00202 DSDPCALLOC2(&ds->an,double,tnnz,&info);DSDPCHKERR(info); 00203 for (i=0;i<tnnz;i++) ds->col[i]=cols[i]; 00204 info=DSDPDSSparseInitializeOpsU(&tdsdsopsu); DSDPCHKERR(info); 00205 *dsmatops=&tdsdsopsu; 00206 *dsmat=(void*)ds; 00207 DSDPFunctionReturn(0); 00208 }