DSDP
dsdp.c
Go to the documentation of this file.
00001 #include "mex.h"
00002 #include "dsdp5.h"
00003 #include <math.h>
00004 
00005 /* #define K_IN prhs[0] */
00006 #define CA_IN prhs[1]
00007 #define B_IN prhs[0]
00008 #define PARS_IN prhs[2]
00009 #define Y_IN prhs[3]
00010 #define STAT_OUT plhs[0]
00011 #define Y_OUT plhs[1]
00012 #define X_OUT plhs[2]
00013 static int DSDPPrintStats2(DSDP, void*);
00014 static int CountNonzeroMatrices(double*,int,int,int,int*,int*,int*);
00015 static int CheckForConstantMat(double*,int,int);
00016 
00017 static int printlevel=10;
00018 #define CHKERR(a)  { if (a){mexWarnMsgTxt("DSDP Numerical Error"); } }
00019 
00020 #define mexErrMsgTxt(a); mexPrintf("Error: "); mexWarnMsgTxt(a); return;
00021 
00035 void mexFunction(int nlhs, mxArray *plhs[],
00036                  int nrhs, const mxArray *prhs[]){
00037   
00038   mxArray  *CA_cell_pr,*X_cell_pr;
00039   const mxArray  *OPTIONS_FIELD;
00040   mxArray  *STAT_FIELD;
00041   int     i,j,k,ii,itmp,index,info;
00042   int    *air, *ajc, *air2, *ajc2, *str,*iptr;
00043   int   nvars,nb,mC,nC,m1,n1,m2,n2;
00044   int   nsubs=2, subs[2];
00045   int   iscellCA;
00046   int   its,reuse=4,print_info=1,printsummary=0;
00047   int   ijnnz,spot,n,nn=0,nzmats,vecn;
00048   int it1,it2;
00049   int nsdpblocks=0,sdpblockj=0,sdpnmax=1,lpnmax=1,stat1=1,xmaker=0;
00050   int sspot,nsubblocks,blockj;
00051   int jj,tnnz,tnnz2;
00052   int maxit=1000,fastblas=1,rpos=0,drho=1,iloginfo=0,aggressive=0;
00053   double penalty=1e8,rho=4,zbar=1e10,cc=0,r0=-1,mu0=-1,ylow,yhigh,gaptol=1e-6,pnormtol=1e30;
00054   double maxtrust=1e30,steptol=0.01,inftol=1e-8,lpb=1.0,dbound=1e20,infptol=1e-4;
00055   double dtmp,pstep,dstep,pnorm,mu;
00056   double *blockn,datanorm[3];
00057   double *aval,*aval2,*bval,*yout,*y0=0,*xout,*stat;
00058   double pobj,dobj,dinf;
00059   DSDP dsdp;
00060   SDPCone  sdpcone=0;
00061   DSDPTerminationReason reason;
00062   DSDPSolutionType pdfeasible;
00063   LPCone   lpcone=0;
00064   BCone bcone=0;
00065   char conetype[30];
00066   int nfields=25;
00067   const char *fnames[25]={"stype","obj","pobj","dobj","stopcode","termcode","iterates","r","mu",
00068                           "pstep","dstep","pnorm","gaphist","infeashist","errors",
00069                           "datanorm","ynorm","boundy","penalty","tracex","reuse","rho","xy","xdy","xmu"};
00070   
00071   if (nrhs < 2){
00072     mexErrMsgTxt("Two input arguments required.  See help for details. ");}
00073   if (nrhs > 4){
00074     mexErrMsgTxt("Fewer input arguments required.  See help for details. ");}
00075   if (nlhs < 2){
00076     mexErrMsgTxt("Two output arguments required.  See help for details. ");}
00077   if (nlhs > 3){
00078     mexErrMsgTxt("Fewer output arguments required.  See help for details. ");}
00079 
00080   if (!mxIsDouble(B_IN) || mxIsSparse(B_IN)){
00081     mexErrMsgTxt("DSDP: 1ST input must be a dense vector of doubles"); }
00082   nvars = mxGetM(B_IN);
00083   nb = mxGetN(B_IN);    
00084   if (nb > 1){
00085     mexErrMsgTxt("DESP: 1ST input must be a column vector"); }
00086   
00087   iscellCA = mxIsCell(CA_IN); 
00088   if (!iscellCA){
00089     mexErrMsgTxt("DSDP: 2ND input must be a cell array"); }
00090   mC = mxGetM(CA_IN);
00091   nC = mxGetN(CA_IN); 
00092   if (nC != 3){
00093     mexErrMsgTxt("DSDP: dimension of 2ND cell array is p x 3");}
00094 
00095   if (nrhs >2){
00096     if(!mxIsStruct(PARS_IN)){
00097       mexErrMsgTxt("3RD input `OPTIONS' should be a structure.");
00098     }
00099   }
00100 
00101   if (nrhs>3){
00102     if (!mxIsDouble(Y_IN) || mxIsSparse(Y_IN)){
00103       mexErrMsgTxt("DSDP: 4TH input must be a dense vector of doubles"); }
00104     m1 = mxGetM(Y_IN);
00105     n1 = mxGetN(Y_IN);
00106     if (m1 != nvars || n1 != nb){
00107       mexErrMsgTxt("DSDP: dimensions of 1ST and 4TH input not compatible");}
00108     y0=mxGetPr(Y_IN);
00109     if (!y0){
00110        mexErrMsgTxt("DSDP: Cannot read 4TH argument");}
00111   } else y0=0;
00112   
00113   
00114   /* Check data */
00115   for (j=0; j<mC; j++){
00116     subs[0] = j; subs[1] = 0;
00117     index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00118     CA_cell_pr = mxGetCell(CA_IN,index); 
00119     if (!CA_cell_pr){
00120       mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,1);
00121       mexErrMsgTxt("DSDP: Empty Cell. Missing String"); }
00122     if (!mxIsChar(CA_cell_pr)){
00123       mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00124       mexErrMsgTxt("DSDP: First column of cells in 2ND argument are a string to determine cone type."); }
00125     mxGetString(CA_cell_pr,conetype,20); 
00126     
00127     subs[0] = j; subs[1] = 1;
00128     index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00129     CA_cell_pr = mxGetCell(CA_IN,index); 
00130     if (!CA_cell_pr){
00131       mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,2);
00132       mexErrMsgTxt("DSDP: Empty Cell. Provide dimension of block."); }
00133     if (mxIsSparse(CA_cell_pr)){
00134       mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00135       mexErrMsgTxt("DSDP: Second column in 2ND argument must be a dense array of scalars that specify dimension."); }
00136     if (!mxIsDouble(CA_cell_pr)){
00137       mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00138       mexErrMsgTxt("DSDP: Second column in 2ND argument must specify dimension."); }
00139     aval=mxGetPr(CA_cell_pr);
00140 
00141     subs[0] = j; subs[1] = 2;
00142     index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00143     CA_cell_pr = mxGetCell(CA_IN,index); 
00144     if (!CA_cell_pr){
00145       mexPrintf("??? DSDP DATA ERROR: Cell: %d, Column: %d\n",j+1,3);
00146       mexErrMsgTxt("DSDP: Empty Cell. Provide sparse data matrix."); }
00147     if (!mxIsSparse(CA_cell_pr) || !mxIsDouble(CA_cell_pr)){ 
00148       mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00149       mexErrMsgTxt("DSDP: Third column in 2ND argument must be a real sparse data matrix."); }
00150     
00151     if (strcmp(conetype,"SDP")==0){ 
00152       subs[0] = j; subs[1] = 1;
00153       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00154       CA_cell_pr = mxGetCell(CA_IN,index); 
00155       aval=mxGetPr(CA_cell_pr);
00156       it1 = mxGetM(CA_cell_pr);
00157       it2 = mxGetN(CA_cell_pr); 
00158       if (it1!=1 && it2!=1){
00159         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00160         mexErrMsgTxt("DSDP: Use a dense row vector in the second column in 2ND of the argument.");}
00161 
00162       subs[0] = j; subs[1] = 2;
00163       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00164       CA_cell_pr = mxGetCell(CA_IN,index); 
00165       m1 = mxGetN(CA_cell_pr)-1;
00166       if (m1 != nvars){
00167         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00168         mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have number of columns equal to number of variables+1.");}
00169       vecn = mxGetM(CA_cell_pr); 
00170       for (tnnz=0,i=0;i<it1*it2;i++){n=(int)aval[i]; tnnz=tnnz+n*(n+1)/2;}
00171       if ( tnnz != vecn){
00172         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00173         mexErrMsgTxt("DSDP: Check Dimensions:  The columns of A and C cannot be converted into square matrices");}
00174       nsdpblocks=nsdpblocks+it1*it2;
00175     } else if (strcmp(conetype,"LP")==0){ 
00176 
00177       subs[0] = j; subs[1] = 1;
00178       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00179       CA_cell_pr = mxGetCell(CA_IN,index); 
00180       it1 = mxGetM(CA_cell_pr);
00181       it2 = mxGetN(CA_cell_pr); 
00182       /* THe sum of the dimensions should equal the number of constraints */
00183       
00184       subs[0] = j; subs[1] = 2;
00185       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00186       CA_cell_pr = mxGetCell(CA_IN,index); 
00187       if (!mxIsSparse(CA_cell_pr)){ 
00188         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00189         mexErrMsgTxt("DSDP: Matrices in the third column in 2ND argument must be sparse."); }     
00190       m1 = mxGetN(CA_cell_pr)-1;
00191       if (m1 != nvars){
00192         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00193         mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have number of columns equal to number of variables+1.");}
00194     } else if (strcmp(conetype,"LB")==0 || strcmp(conetype,"UB")==0 ){ 
00195       
00196       subs[0] = j; subs[1] = 2;
00197       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00198       CA_cell_pr = mxGetCell(CA_IN,index); 
00199       if (mxIsSparse(CA_cell_pr)){ 
00200         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00201         mexErrMsgTxt("DSDP: Row vector in the third column in 2ND argument must be full."); }     
00202       it1 = mxGetM(CA_cell_pr);
00203       it2 = mxGetN(CA_cell_pr); 
00204       if (it2 != 1){
00205         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00206         mexErrMsgTxt("DSDP: The matrix in the third column in 2ND argument must have a single column of bounds.");}
00207       if (it1 != nvars){
00208         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00209         mexErrMsgTxt("DSDP: The column matrix in the third column in 2ND argument must contain a bound for each y variable.");}
00210 
00211     } else if (strcmp(conetype,"FIXED")==0){ 
00212       int dim1;
00213       subs[0] = j; subs[1] = 1;
00214       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00215       CA_cell_pr = mxGetCell(CA_IN,index); 
00216       if (mxIsSparse(CA_cell_pr)){ 
00217         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00218         mexErrMsgTxt("DSDP: Vector in the third column in 2ND argument must be full."); }     
00219       it1 = mxGetM(CA_cell_pr);
00220       it2 = mxGetN(CA_cell_pr); 
00221       dim1=it1*it2;
00222       if (it1 != 1){
00223         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00224         mexErrMsgTxt("DSDP: Third column in 2ND argument must have 1 row,");}
00225       subs[0] = j; subs[1] = 2;
00226       index = mxCalcSingleSubscript(CA_IN,nsubs,subs);
00227       CA_cell_pr = mxGetCell(CA_IN,index); 
00228       if (mxIsSparse(CA_cell_pr)){ 
00229         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00230         mexErrMsgTxt("DSDP: Vector in the third column in 2ND argument must be full."); }
00231       it1 = mxGetM(CA_cell_pr);
00232       it2 = mxGetN(CA_cell_pr); 
00233       if (it1 != 1){
00234         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00235         mexErrMsgTxt("DSDP: Third column in 2ND argument must have 1 row,");}
00236       if (it2 != dim1){
00237         mexPrintf("??? DSDP DATA ERROR: Cell: %d \n",j+1);
00238         mexErrMsgTxt("DSDP: Secord and third column must have same dimension,");}
00239     } else {
00240       mexPrintf("??? DSDP DATA ERROR: Cell: %d, Conetype: %s \n",j+1,conetype);
00241       mexErrMsgTxt("DSDP: Unknown Cone type in 2ND argument. Try 'SDP' or 'LP' or 'Bounds'. ");
00242     }
00243   }
00244   
00245   
00246   /* Create output arrays */
00247   if (nlhs>2){
00248     if (X_OUT != NULL) mxDestroyArray(X_OUT) ;
00249     X_OUT = mxCreateCellMatrix(mC,1);
00250   }
00251   if (Y_OUT != NULL) mxDestroyArray(Y_OUT) ;
00252   Y_OUT = mxCreateDoubleMatrix(nvars, 1, mxREAL) ;
00253   
00254   /* Create the Solver */
00255   info = DSDPCreate(nvars,&dsdp); CHKERR(info);
00256   info = DSDPCreateSDPCone(dsdp,nsdpblocks,&sdpcone); CHKERR(info);
00257 
00258   /* Set Dual Objective Vector */
00259   bval=mxGetPr(B_IN);
00260   if (!bval){ mexErrMsgTxt("DSDP: Problems with 1ST argument");}
00261   for (i=0;i<nvars;i++){info=DSDPSetDualObjective(dsdp,i+1,bval[i]); CHKERR(info);}
00262   
00263   /* Set Matrix Data */
00264   for (j=0; j<mC; j++){  /* Begin Block */
00265     subs[0] = j; subs[1] = 0;
00266     index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00267     CA_cell_pr = mxGetCell(CA_IN,index); 
00268     mxGetString(CA_cell_pr,conetype,20); 
00269     if (strcmp(conetype,"SDP")==0){ 
00270 
00271       subs[0] = j; subs[1] = 1;
00272       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00273       CA_cell_pr = mxGetCell(CA_IN,index); 
00274       it1 = mxGetM(CA_cell_pr);
00275       it2 = mxGetN(CA_cell_pr); 
00276       blockn=mxGetPr(CA_cell_pr);
00277       nsubblocks=it1*it2;
00278 
00279       subs[0] = j; subs[1] = 2;
00280       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00281       CA_cell_pr = mxGetCell(CA_IN,index); 
00282       aval=mxGetPr(CA_cell_pr); air =mxGetIr(CA_cell_pr); ajc =mxGetJc(CA_cell_pr);
00283       if (!aval||!air||!ajc)
00284         { mexErrMsgTxt("DSDP: Problems with 2ND argument");}
00285       for (tnnz=0,jj=0;jj<nsubblocks;jj++){n=(int)blockn[jj];tnnz+=n*(n+1)/2;}
00286       if (nlhs>2){
00287         subs[0] = j; subs[1] = 0;
00288         index = mxCalcSingleSubscript(X_OUT,nsubs,subs);
00289         X_cell_pr = mxCreateDoubleMatrix(tnnz,1,mxREAL);
00290         mxSetCell(X_OUT,index,X_cell_pr);
00291         xout=mxGetPr(X_cell_pr);
00292         if (tnnz>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
00293       }
00294       
00295       for (ii=0; ii<=nvars; ii++){ /* Begin Variable matrix constraints */
00296         i=ii+1;
00297         if (i==nvars+1) i=0;
00298         tnnz=0; spot=ajc[ii]; blockj=sdpblockj;
00299         for (jj=0;jj<nsubblocks;jj++){
00300           n=(int)blockn[jj];
00301           if (sdpnmax<n) sdpnmax=n;
00302           if (ii==0){
00303             nn+=n;
00304             info=SDPConeSetBlockSize(sdpcone,blockj,n); CHKERR(info);
00305             info=SDPConeUsePackedFormat(sdpcone,blockj); CHKERR(info);
00306 
00307             if (nlhs>2){info=SDPConeSetXArray(sdpcone,blockj,n,xout+tnnz,(n*n+n)/2); CHKERR(info);} 
00308             info=CountNonzeroMatrices(blockn,nsubblocks,jj,nvars, air, ajc, &nzmats); CHKERR(info);
00309             info=SDPConeSetSparsity(sdpcone,blockj,nzmats); CHKERR(info);
00310             if (stat1<nzmats)stat1=nzmats;
00311           }
00312           for (tnnz2=tnnz+n*(n+1)/2,ijnnz=0;ijnnz<ajc[ii+1]-spot && air[spot+ijnnz]<tnnz2;ijnnz++){}
00313           if ( ijnnz==0 ){     /*  info=DSDPSetZeroMat(dsdp,sdpblockj,i,n); */
00314           } else if(ijnnz==n*(n+1)/2){ /* check for dense matrix  */
00315             if (CheckForConstantMat(aval+spot,ijnnz,n)){
00316               info=SDPConeSetConstantMat(sdpcone,blockj,i,n,aval[spot]); CHKERR(info);
00317             } else {
00318               info=SDPConeSetADenseVecMat(sdpcone,blockj,i,n,1.0,aval+spot,ijnnz); CHKERR(info);
00319             }
00320           } else {
00321             info=SDPConeSetASparseVecMat(sdpcone,blockj,i,n,1.0,tnnz,air+spot,aval+spot,ijnnz); CHKERR(info);
00322           }
00323           tnnz+=n*(n+1)/2; spot+=ijnnz; blockj++;
00324         }
00325       }  /* End Matrices in block */
00326       sdpblockj=sdpblockj+nsubblocks;
00327 
00328     } else if (strcmp(conetype,"LB")==0 || strcmp(conetype,"UB")==0 ){
00329       subs[0] = j; subs[1] = 2;
00330       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00331       CA_cell_pr = mxGetCell(CA_IN,index);
00332       aval=mxGetPr(CA_cell_pr);
00333       info=DSDPCreateBCone(dsdp,&bcone);  CHKERR(info);
00334       info=BConeAllocateBounds(bcone,nvars); CHKERR(info);
00335       for (i=0;i<nvars;i++){
00336         if (strcmp(conetype,"LB")==0){
00337           info=BConeSetLowerBound(bcone,i+1,aval[i]); CHKERR(info);
00338         } else {
00339           info=BConeSetUpperBound(bcone,i+1,aval[i]); CHKERR(info);
00340         } 
00341       }
00342       if (nlhs>2){
00343         subs[0] = j; subs[1] = 0;
00344         index = mxCalcSingleSubscript(X_OUT,nsubs,subs); 
00345         X_cell_pr = mxCreateDoubleMatrix(nvars,1,mxREAL);
00346         mxSetCell(X_OUT,index,X_cell_pr);
00347         aval2=mxGetPr(X_cell_pr);
00348         info=BConeSetXArray(bcone,aval2,nvars); CHKERR(info);
00349       }
00350     } else if (strcmp(conetype,"LP")==0){
00351       subs[0] = j; subs[1] = 2;
00352       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00353       CA_cell_pr = mxGetCell(CA_IN,index); 
00354       n = mxGetM(CA_cell_pr); 
00355       if (lpnmax<n) lpnmax=n;
00356       nn+=n;
00357       aval=mxGetPr(CA_cell_pr); air =mxGetIr(CA_cell_pr); ajc =mxGetJc(CA_cell_pr);
00358       if (!aval||!air||!ajc)
00359         { mexErrMsgTxt("DSDP: Problems with 2ND argument");}
00360       info=DSDPCreateLPCone(dsdp,&lpcone); CHKERR(info);
00361       info=LPConeSetData2(lpcone,n,ajc,air,aval); CHKERR(info);
00362       if (nlhs>2){
00363         subs[0] = j; subs[1] = 0;
00364         index = mxCalcSingleSubscript(X_OUT,nsubs,subs); 
00365         X_cell_pr = mxCreateDoubleMatrix(n,1,mxREAL);
00366         mxSetCell(X_OUT,index,X_cell_pr);
00367         xout=mxGetPr(X_cell_pr);
00368         if (n>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
00369         info=LPConeSetXVec(lpcone,xout,n); CHKERR(info);
00370       }
00371     
00372     } else if (strcmp(conetype,"FIXED")==0){
00373       int vari;
00374       subs[0] = j; subs[1] = 1;
00375       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00376       CA_cell_pr = mxGetCell(CA_IN,index);
00377       aval=mxGetPr(CA_cell_pr);
00378       it1 = mxGetM(CA_cell_pr);
00379       it2 = mxGetN(CA_cell_pr); 
00380       subs[0] = j; subs[1] = 2;
00381       index = mxCalcSingleSubscript(CA_IN,nsubs,subs); 
00382       CA_cell_pr = mxGetCell(CA_IN,index);
00383       aval2=mxGetPr(CA_cell_pr);
00384       if (nlhs>2){
00385         subs[0] = j; subs[1] = 0;
00386         index = mxCalcSingleSubscript(X_OUT,nsubs,subs); 
00387         X_cell_pr = mxCreateDoubleMatrix(it1*it2,1,mxREAL);
00388         mxSetCell(X_OUT,index,X_cell_pr);
00389         xout=mxGetPr(X_cell_pr);
00390         if (it1*it2>0 && !xout){ mexErrMsgTxt("DSDP: Cannot create array. Out of Memory");}
00391       }  else {xout=0;}
00392       info=DSDPSetFixedVariables(dsdp,aval,aval2,xout,it1*it2); CHKERR(info);
00393       for (i=0;i<it1*it2;i++){
00394         /*
00395         vari=(int)aval[i];
00396         printf("FixedVariable %d to %4.4e\n",vari,aval2[i]);
00397         info=DSDPSetFixedVariable(dsdp,vari,aval2[i]); CHKERR(info); 
00398         */
00399       }
00400     }    
00401   } /* End Block */
00402 
00403   /* Set initial point */
00404   if (y0){
00405     for (i=0;i<nvars;i++){ info = DSDPSetY0(dsdp,i+1,y0[i]); CHKERR(info);}
00406   } 
00407 
00408   reuse=(nvars-2)/sdpnmax; 
00409   if (nvars<50 && reuse==0) reuse=1;
00410   if (reuse>=1) reuse++;
00411   reuse=reuse*reuse;
00412   if (nvars<2000 && reuse>10) reuse=10;
00413   if (reuse>12) reuse=12;
00414   info=DSDPReuseMatrix(dsdp,reuse); CHKERR(info);
00415 
00416   info=DSDPGetPPObjective(dsdp,&zbar); CHKERR(info);
00417   info=DSDPGetR(dsdp,&r0); CHKERR(info);
00418   info=DSDPGetMaxIts(dsdp,&maxit);  CHKERR(info);
00419   info=DSDPGetPenaltyParameter(dsdp,&penalty);  CHKERR(info);
00420   info=DSDPGetPotentialParameter(dsdp,&rho);  CHKERR(info);
00421   info=DSDPGetDualBound(dsdp,&dbound);  CHKERR(info);
00422   info=DSDPGetGapTolerance(dsdp,&gaptol); CHKERR(info);
00423   info=DSDPGetRTolerance(dsdp,&inftol);  CHKERR(info);
00424   info=DSDPGetBarrierParameter(dsdp,&mu0); CHKERR(info);
00425   info=DSDPGetMaxTrustRadius(dsdp,&maxtrust);  CHKERR(info);
00426   info=DSDPGetStepTolerance(dsdp,&steptol);  CHKERR(info);
00427   info=DSDPGetPTolerance(dsdp,&infptol);  CHKERR(info);
00428   info=DSDPGetDataNorms(dsdp, datanorm); CHKERR(info);
00429   info=DSDPGetYBounds(dsdp,&ylow,&yhigh); CHKERR(info);
00430   if (datanorm[0]==0){info=DSDPSetYBounds(dsdp,-1.0,1.0);CHKERR(info);}
00431   if (nrhs >2){
00432     if(!mxIsStruct(PARS_IN)){
00433       mexErrMsgTxt("Fifth Parameter `OPTIONS' should be a structure.");}
00434 
00435       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"maxit") ){
00436         maxit=(int) mxGetScalar(OPTIONS_FIELD);
00437         info=DSDPSetMaxIts(dsdp,maxit); CHKERR(info);}
00438       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"fastblas") ){
00439         fastblas= (int) mxGetScalar(OPTIONS_FIELD);}
00440       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"print") ){
00441         print_info= (int) mxGetScalar(OPTIONS_FIELD);
00442         printlevel=print_info;
00443       }
00444       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"bigM") ){
00445         rpos=(int)mxGetScalar(OPTIONS_FIELD);
00446         info=DSDPUsePenalty(dsdp,rpos); CHKERR(info);}
00447       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"penalty") ){
00448         penalty=mxGetScalar(OPTIONS_FIELD);
00449         info=DSDPSetPenaltyParameter(dsdp,penalty); CHKERR(info);}
00450       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"rho") ){
00451         rho= mxGetScalar(OPTIONS_FIELD);
00452         info=DSDPSetPotentialParameter(dsdp,rho); CHKERR(info);}
00453       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dynamicrho") ){
00454         drho= (int)mxGetScalar(OPTIONS_FIELD);
00455         info=DSDPUseDynamicRho(dsdp,drho); CHKERR(info);}
00456       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"zbar") ){
00457         zbar= mxGetScalar(OPTIONS_FIELD);
00458         info=DSDPSetZBar(dsdp,zbar); CHKERR(info);}
00459       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dual_bound") ){
00460         dbound= mxGetScalar(OPTIONS_FIELD);
00461         info=DSDPSetDualBound(dsdp,dbound); CHKERR(info);}
00462       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"reuse") ){
00463         reuse= (int)mxGetScalar(OPTIONS_FIELD);
00464         info=DSDPReuseMatrix(dsdp,reuse); CHKERR(info);}
00465       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"gaptol") ){
00466         gaptol= mxGetScalar(OPTIONS_FIELD);
00467         info=DSDPSetGapTolerance(dsdp,gaptol);CHKERR(info);}
00468       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"lp_barrier") ){
00469         lpb= mxGetScalar(OPTIONS_FIELD);
00470         if (lpb<0.1) lpb=0.1;
00471         if (lpcone){info=LPConeScaleBarrier(lpcone,lpb); CHKERR(info);}}
00472       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"lpb") ){
00473         lpb= mxGetScalar(OPTIONS_FIELD);
00474         if (lpb<0.1) lpb=0.1;
00475         if (lpcone){info=LPConeScaleBarrier(lpcone,lpb); CHKERR(info);}}
00476       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"cc") ){
00477         cc= mxGetScalar(OPTIONS_FIELD);
00478         info=DSDPAddObjectiveConstant(dsdp,cc); CHKERR(info);}
00479       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"inftol") ){
00480         inftol= mxGetScalar(OPTIONS_FIELD);
00481         info=DSDPSetRTolerance(dsdp,inftol); CHKERR(info);}
00482       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"infptol") ){
00483         infptol= mxGetScalar(OPTIONS_FIELD);
00484         info=DSDPSetPTolerance(dsdp,infptol);CHKERR(info);}
00485       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"pnormtol") ){
00486         pnormtol= mxGetScalar(OPTIONS_FIELD);
00487         info=DSDPSetPNormTolerance(dsdp,pnormtol);CHKERR(info);}
00488       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"boundy") ){
00489         yhigh= fabs(mxGetScalar(OPTIONS_FIELD)); ylow=-yhigh;
00490         info=DSDPSetYBounds(dsdp,ylow,yhigh);CHKERR(info);}
00491       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"r0") ){
00492         r0= mxGetScalar(OPTIONS_FIELD); 
00493         info=DSDPSetR0(dsdp,r0); CHKERR(info);}
00494       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"mu0") ){
00495         mu0=mxGetScalar(OPTIONS_FIELD);
00496         info=DSDPSetBarrierParameter(dsdp,mu0);CHKERR(info);}
00497       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"maxtrustradius")){
00498         maxtrust= mxGetScalar(OPTIONS_FIELD);
00499         info=DSDPSetMaxTrustRadius(dsdp,maxtrust); CHKERR(info);}
00500       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"steptol") ){
00501         steptol = mxGetScalar(OPTIONS_FIELD);
00502         info=DSDPSetStepTolerance(dsdp,steptol); CHKERR(info);}
00503       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dobjmin") ){
00504         dtmp= mxGetScalar(OPTIONS_FIELD);
00505         info = DSDPSetDualLowerBound(dsdp,dtmp);}
00506       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"dloginfo") ){
00507         iloginfo= (int) mxGetScalar(OPTIONS_FIELD);
00508         info=DSDPLogInfoAllow(iloginfo,0);}
00509       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"logtime") ){
00510         printsummary= (int) mxGetScalar(OPTIONS_FIELD);}
00511       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"printproblem") ){
00512         info=DSDPPrintData(dsdp,sdpcone,lpcone);}
00513       if ( OPTIONS_FIELD = mxGetField(PARS_IN,0,"xmaker") ){
00514         xmaker = (int) mxGetScalar(OPTIONS_FIELD);}
00515       /*
00516         if( OPTIONS_FIELD = mxGetField(PARS_IN,0,"maxlanczos") ){
00517         itmp= (int) mxGetScalar(OPTIONS_FIELD);
00518         info=DSDPSetLanczosIterations(dsdp,itmp);}
00519       */
00520       
00521   }
00522 
00523   info = DSDPSetMonitor(dsdp,DSDPPrintStats2,0); CHKERR(info);
00524 
00525   info = DSDPSetup(dsdp);
00526   if (info){
00527     mexErrMsgTxt("DSDP: Setup Error, Probably out of memory");}
00528 
00529 
00530   info = DSDPSolve(dsdp); CHKERR(info);
00531   info = DSDPStopReason(dsdp,&reason); 
00532   if (reason!=DSDP_INFEASIBLE_START){
00533     info=DSDPComputeX(dsdp);CHKERR(info);
00534   }
00535   info = DSDPGetSolutionType(dsdp,&pdfeasible);  CHKERR(info);
00536 
00537   if (printsummary){ DSDPEventLogSummary();}
00538 
00539   if (info){
00540     mexErrMsgTxt("DSDP: Numerical error");}
00541 
00542   if ( reason == DSDP_INFEASIBLE_START){
00543     mexErrMsgTxt("DSDP Terminated Due to Infeasible Starting Point\n");
00544   } else if (print_info){
00545 
00546     if (reason == DSDP_CONVERGED)
00547       mexPrintf("DSDP Converged. \n");
00548     else if ( reason == DSDP_UPPERBOUND )
00549       mexPrintf("DSDP Converged: Dual Objective exceeds its bound\n");
00550     else if ( reason == DSDP_SMALL_STEPS )
00551       mexPrintf("DSDP Terminated Due to Small Steps\n");
00552     else if ( reason == DSDP_MAX_IT)
00553       mexPrintf("DSDP Terminated Due Maximum Number of Iterations\n");
00554     else if ( reason == DSDP_USER_TERMINATION)
00555       mexPrintf("DSDP Terminated By User\n");
00556     else if ( reason == DSDP_INFEASIBLE_START)
00557       mexPrintf("DSDP Terminated Due to Infeasible Starting Point\n");
00558     else 
00559       mexPrintf("DSDP Finished.\n");
00560   }
00561 
00562   if (pdfeasible==DSDP_UNBOUNDED){
00563     mexPrintf("DSDP: Dual Unbounded, Primal Infeasible\n");
00564   } else if (pdfeasible==DSDP_INFEASIBLE){
00565     mexPrintf("DSDP: Primal Unbounded, Dual Infeasible\n");
00566   }
00567 
00568   /* Set the dual solution */  
00569   yout=mxGetPr(Y_OUT);
00570   info = DSDPGetY(dsdp,yout,nvars); CHKERR(info);
00571   if (info){
00572     mexErrMsgTxt("DSDP: Numerical error");}
00573   
00574   
00575   /* Output statistics */
00576   if (STAT_OUT != NULL) mxDestroyArray(STAT_OUT) ;
00577   subs[0] = 1;  subs[1] = 1;
00578   STAT_OUT = mxCreateStructArray(2,subs,nfields,fnames);
00579   info= DSDPGetDObjective(dsdp,&dobj);  CHKERR(info);
00580   info= DSDPGetPObjective(dsdp,&pobj);  CHKERR(info);
00581   info= DSDPGetR(dsdp,&dinf);  CHKERR(info);
00582   info= DSDPStopReason(dsdp,&reason);  CHKERR(info);
00583   info= DSDPGetIts(dsdp,&its);  CHKERR(info);
00584   info= DSDPGetBarrierParameter(dsdp,&mu); CHKERR(info);
00585   info= DSDPGetStepLengths(dsdp,&pstep,&dstep); CHKERR(info);
00586   info= DSDPGetPnorm(dsdp,&pnorm); CHKERR(info);
00587   info= DSDPGetBarrierParameter(dsdp,&mu); CHKERR(info);
00588   info= DSDPGetYBounds(dsdp,&ylow,&yhigh); CHKERR(info);
00589 
00590   if (pdfeasible==DSDP_UNBOUNDED){
00591     STAT_FIELD = mxCreateString("Unbounded");
00592     mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
00593   } else if (pdfeasible==DSDP_INFEASIBLE){
00594     STAT_FIELD = mxCreateString("Infeasible");
00595     mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
00596   } else {
00597     STAT_FIELD = mxCreateString("PDFeasible");
00598     mxSetField(STAT_OUT,0,fnames[0],STAT_FIELD);
00599   }
00600 
00601   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00602   stat=mxGetPr(STAT_FIELD); stat[0]=dobj;
00603   mxSetField(STAT_OUT,0,fnames[1],STAT_FIELD);
00604 
00605   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00606   stat=mxGetPr(STAT_FIELD); stat[0]=pobj;
00607   mxSetField(STAT_OUT,0,fnames[2],STAT_FIELD);
00608   
00609   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00610   stat=mxGetPr(STAT_FIELD); stat[0]=dobj;
00611   mxSetField(STAT_OUT,0,fnames[3],STAT_FIELD);
00612 
00613   if (reason==DSDP_CONVERGED){
00614     STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00615     stat=mxGetPr(STAT_FIELD); stat[0]=0;
00616     mxSetField(STAT_OUT,0,fnames[4],STAT_FIELD);
00617   } else {
00618     STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00619     stat=mxGetPr(STAT_FIELD); stat[0]=(double)reason;
00620     if (stat[0]==0) stat[0]=-1;
00621     mxSetField(STAT_OUT,0,fnames[4],STAT_FIELD);
00622   }
00623   
00624   if (reason==DSDP_CONVERGED){ /* For YALMIP */
00625     STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00626     stat=mxGetPr(STAT_FIELD); stat[0]=0;
00627     if (pdfeasible==DSDP_UNBOUNDED){ stat[0]=1;} 
00628     if (pdfeasible==DSDP_INFEASIBLE){ stat[0]=2;}
00629     mxSetField(STAT_OUT,0,fnames[5],STAT_FIELD);
00630   } else {
00631     STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00632     stat=mxGetPr(STAT_FIELD); stat[0]=-1;
00633     if (reason==DSDP_INFEASIBLE_START)stat[0]=-6;
00634     if (reason==DSDP_UPPERBOUND)stat[0]=-5;
00635     if (reason==DSDP_MAX_IT)stat[0]=-3;
00636     mxSetField(STAT_OUT,0,fnames[5],STAT_FIELD);
00637   }
00638 
00639   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00640   stat=mxGetPr(STAT_FIELD); stat[0]=(double) its;
00641   mxSetField(STAT_OUT,0,fnames[6],STAT_FIELD);
00642   
00643   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00644   stat=mxGetPr(STAT_FIELD); stat[0]=dinf;
00645   mxSetField(STAT_OUT,0,fnames[7],STAT_FIELD);
00646     
00647   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00648   stat=mxGetPr(STAT_FIELD); stat[0]=mu;
00649   mxSetField(STAT_OUT,0,fnames[8],STAT_FIELD);
00650   
00651   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00652   stat=mxGetPr(STAT_FIELD); stat[0]=pstep;
00653   mxSetField(STAT_OUT,0,fnames[9],STAT_FIELD);
00654   
00655   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00656   stat=mxGetPr(STAT_FIELD); stat[0]=dstep;
00657   mxSetField(STAT_OUT,0,fnames[10],STAT_FIELD);
00658   
00659   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00660   stat=mxGetPr(STAT_FIELD); stat[0]=pnorm;
00661   mxSetField(STAT_OUT,0,fnames[11],STAT_FIELD);
00662   
00663   itmp=100; if (its < itmp) itmp=its;
00664   STAT_FIELD = mxCreateDoubleMatrix(itmp, 1, mxREAL) ;
00665   stat=mxGetPr(STAT_FIELD); 
00666   info= DSDPGetGapHistory(dsdp,stat,itmp);  CHKERR(info);
00667   mxSetField(STAT_OUT,0,fnames[12],STAT_FIELD);
00668   
00669   STAT_FIELD = mxCreateDoubleMatrix(itmp, 1, mxREAL) ;
00670   stat=mxGetPr(STAT_FIELD); 
00671   info= DSDPGetRHistory(dsdp,stat,itmp);  CHKERR(info);
00672   mxSetField(STAT_OUT,0,fnames[13],STAT_FIELD);
00673 
00674   STAT_FIELD = mxCreateDoubleMatrix(6, 1, mxREAL) ;
00675   stat=mxGetPr(STAT_FIELD); 
00676   info = DSDPGetFinalErrors(dsdp,stat);  CHKERR(info);
00677   mxSetField(STAT_OUT,0,fnames[14],STAT_FIELD);
00678 
00679   STAT_FIELD = mxCreateDoubleMatrix(3, 1, mxREAL) ;
00680   stat=mxGetPr(STAT_FIELD); 
00681   info = DSDPGetDataNorms(dsdp,stat);  CHKERR(info);
00682   dtmp=stat[0];stat[0]=stat[1];stat[1]=dtmp;
00683   mxSetField(STAT_OUT,0,fnames[15],STAT_FIELD);
00684 
00685   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00686   stat=mxGetPr(STAT_FIELD); 
00687   info = DSDPGetYMaxNorm(dsdp,stat);  CHKERR(info);
00688   mxSetField(STAT_OUT,0,fnames[16],STAT_FIELD);
00689 
00690   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00691   stat=mxGetPr(STAT_FIELD); 
00692   stat[0]=yhigh;
00693   mxSetField(STAT_OUT,0,fnames[17],STAT_FIELD);
00694 
00695   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00696   stat=mxGetPr(STAT_FIELD); 
00697   info = DSDPGetPenaltyParameter(dsdp,stat);  CHKERR(info);
00698   mxSetField(STAT_OUT,0,fnames[18],STAT_FIELD);
00699 
00700   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00701   stat=mxGetPr(STAT_FIELD); 
00702   info = DSDPGetTraceX(dsdp,stat);  CHKERR(info);
00703   mxSetField(STAT_OUT,0,fnames[19],STAT_FIELD);
00704 
00705   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00706   info = DSDPGetReuseMatrix(dsdp,&its);  CHKERR(info);
00707   stat=mxGetPr(STAT_FIELD); stat[0]=(double) its;
00708   mxSetField(STAT_OUT,0,fnames[20],STAT_FIELD);
00709   
00710   STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00711   stat=mxGetPr(STAT_FIELD); 
00712   info = DSDPGetPotentialParameter(dsdp,stat);  CHKERR(info);
00713   mxSetField(STAT_OUT,0,fnames[21],STAT_FIELD);
00714 
00715   if (xmaker){
00716     STAT_FIELD = mxCreateDoubleMatrix(nvars+1, 1, mxREAL) ;
00717     stat=mxGetPr(STAT_FIELD); 
00718     info = DSDPGetYMakeX(dsdp,stat,nvars+1);  CHKERR(info);
00719     mxSetField(STAT_OUT,0,fnames[22],STAT_FIELD);
00720 
00721     STAT_FIELD = mxCreateDoubleMatrix(nvars+1, 1, mxREAL) ;
00722     stat=mxGetPr(STAT_FIELD); 
00723     info = DSDPGetDYMakeX(dsdp,stat,nvars+1);  CHKERR(info);
00724     mxSetField(STAT_OUT,0,fnames[23],STAT_FIELD);
00725 
00726     STAT_FIELD = mxCreateDoubleMatrix(1, 1, mxREAL) ;
00727     stat=mxGetPr(STAT_FIELD); 
00728     info = DSDPGetMuMakeX(dsdp,stat);  CHKERR(info);
00729     mxSetField(STAT_OUT,0,fnames[24],STAT_FIELD);
00730   }
00731   /* Free internal data structure */
00732 
00733   info = DSDPDestroy(dsdp); CHKERR(info);
00734 
00735   return;
00736 } /* main */
00737   
00738 
00739 
00740 #undef __FUNCT__
00741 #define __FUNCT__ "CheckForConstantMat"
00742 static int CheckForConstantMat(double*v,int nnz, int n){
00743   int i;double vv;
00744   if (n<=1){ return 0; }
00745   if (nnz!=(n*n+n)/2){ return 0; }
00746   for (vv=v[0],i=1;i<nnz;i++){
00747     if (v[i]!=vv){ return 0;}
00748   }
00749   return 1;
00750 }
00751 
00752 #undef __FUNCT__
00753 #define __FUNCT__ "CountNonzeroMatrices"
00754  static int CountNonzeroMatrices(double *blocksize, int nblocks, int block, int nvars, int *indd, int*nnnz, int *nnzmats){
00755    int i,j,n,nzmats=0;
00756    int marker1=0,marker2=0;
00757   for (i=0;i<block;i++){n=(int)blocksize[i]; marker1=marker1+n*(n+1)/2;}
00758   n=(int)blocksize[block];  
00759   marker2=marker1+n*(n+1)/2;
00760   for (i=0;i<nvars;i++){
00761     j=nnnz[i];
00762     while (indd[j]<marker1 && j< nnnz[i+1]){j++;}
00763     if (j<nnnz[i+1] && indd[j]<marker2){ nzmats++;}
00764   }
00765   *nnzmats=nzmats;
00766   return 0;
00767 }
00768 
00769 /* ---------------------------------------------------------- */
00770  static int DSDPPrintStats2(DSDP dsdp, void* dummy){
00771    
00772   int    iter,info;
00773   double pobj,dobj,pstp=0,dstp,mu,res,pnorm,pinfeas;
00774   DSDPTerminationReason reason;
00775 
00776   if(printlevel<=0) return(0);
00777 
00778   info = DSDPStopReason(dsdp,&reason);
00779   info = DSDPGetIts(dsdp,&iter);
00780 
00781   if( (reason!=CONTINUE_ITERATING) || ((iter % printlevel)==0)){
00782 
00783     info = DSDPGetDDObjective(dsdp,&dobj);
00784     info = DSDPGetPPObjective(dsdp,&pobj);
00785     info = DSDPGetR(dsdp,&res);
00786     info = DSDPGetPInfeasibility(dsdp,&pinfeas);
00787     info = DSDPGetStepLengths(dsdp,&pstp,&dstp);
00788     info = DSDPGetBarrierParameter(dsdp,&mu);
00789     info = DSDPGetPnorm(dsdp,&pnorm);
00790     
00791     
00792     if (iter==0){
00793       mexPrintf("Iter   PP Objective      DD Objective     PInfeas   DInfeas     Nu     StepLength   Pnrm\n")
00794         ;
00795       mexPrintf("----------------------------------------------------------------------------------------\n")
00796         ;
00797     }
00798     mexPrintf("%-3d %16.8e  %16.8e  %9.1e %9.1e %9.1e",iter,pobj,dobj,pinfeas,res,mu);
00799     mexPrintf("  %4.2f  %4.2f",pstp,dstp);
00800     if (pnorm>1.0e3){
00801       mexPrintf("  %1.0e \n",pnorm);
00802     } else {
00803       mexPrintf("  %5.2f \n",pnorm);
00804     }
00805     fflush(NULL);
00806   }
00807   return(0);
00808 }