分析非线性方面的工具箱--matcont3p1

/* BVP_BPC_jacCC.c MEX file corresponding to BVPjac.m Does the evaluation of the jacobian of the BVP calling syntax: result = BVP_BPC_jacCC(lds.func,x,p,T,pars,nc,lds,gds.period,p2) */ #include<math.h> #include<mex.h> #include<matrix.h> #include<stdio.h> void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* Declarations */ /* ------------ */ double *x,*p,*pars,*pars2, *nc; mxArray *thisfield; int ntst, ncol, nphase, *ActiveParams, ncoords, nfreep, bfreep, tps, *BranchParam; double *upoldp, *mesh, *wt, *wp, *pwi, *BPC_phi1, *BPC_phi2, *BPC_psi; double *dt, *wploc, T; double *pr; mwIndex *ir, *jc; int ncol_coord; mxArray *evalrhs[1000], *jacrhsBR[6], *jacrhs[5]; mxArray *evallhs[1], *jaclhsBR[1], *jaclhs[1]; double *xtmp, *xtmpm, *xtmpn; int filled, elementcounter, remm; int i,j,k,l,l2; /* Indexation variables */ int *range1, *range2, *range3, *range4; double *jac, *jacp, *icjac, *sysjac, *sysjacp, *zero, *ptmp; double *frhs, *frhstmp; double *Tcol, *freepcols, *tempmatrix; double tmpperiod; /* Initializations */ /* --------------- */ /* Retrieve parameters. */ x = mxGetPr(prhs[1]); p = mxGetPr(prhs[2]); pars2 = mxGetPr(prhs[4]); nc = mxGetPr(prhs[5]); /* LDS FIELDS */ thisfield = mxGetFieldByNumber(prhs[6],0,10); nphase = *(mxGetPr(thisfield)); /* Size of one point */ thisfield = mxGetFieldByNumber(prhs[6],0,11); nfreep = mxGetNumberOfElements(thisfield);/* number of free parameters */ ActiveParams = mxCalloc(nfreep,sizeof(int)); frhstmp = mxGetPr(thisfield); for (i=0; i<nfreep; i++) *(ActiveParams+i) = (int)(*(frhstmp+i)); thisfield = mxGetFieldByNumber(prhs[6],0,13); ntst = *(mxGetPr(thisfield)); /* Number of mesh intervals */ thisfield = mxGetFieldByNumber(prhs[6],0,14); ncol = *(mxGetPr(thisfield)); /* Number of collocation points */ thisfield = mxGetFieldByNumber(prhs[6],0,17); tps = *(mxGetPr(thisfield)); thisfield = mxGetFieldByNumber(prhs[6],0,18); ncoords = *(mxGetPr(thisfield)); thisfield = mxGetFieldByNumber(prhs[6],0,22); mesh = mxGetPr(thisfield); /* Current mesh coordinates */ thisfield = mxGetFieldByNumber(prhs[6],0,24); dt = mxGetPr(thisfield); /* Interval widths */ thisfield = mxGetFieldByNumber(prhs[6],0,25); upoldp = mxGetPr(thisfield); /* Derivative of cycle at old mesh coordinates */ thisfield = mxGetFieldByNumber(prhs[6],0,29); wt = mxGetPr(thisfield); /* Weights of collocation points */ thisfield = mxGetFieldByNumber(prhs[6],0,31); T = *(mxGetPr(thisfield)); /* period */ thisfield = mxGetFieldByNumber(prhs[6],0,34); ncol_coord = *(mxGetPr(thisfield)); thisfield = mxGetFieldByNumber(prhs[6],0,39); wp = mxGetPr(thisfield); /* Derivative weights of collocation points */ /* Kronecker product of the derivative weights and the identity matrix */ wploc = mxCalloc(mxGetN(thisfield)*mxGetM(thisfield),sizeof(double)); thisfield = mxGetFieldByNumber(prhs[6],0,41); pwi = mxGetPr(thisfield); /* Extension of weights */ thisfield = mxGetFieldByNumber(prhs[6],0,86); bfreep = mxGetNumberOfElements(thisfield);/* number of branch parameters */ BranchParam = mxCalloc(bfreep,sizeof(int)); frhstmp = mxGetPr(thisfield); for (i=0; i<bfreep; i++) *(BranchParam+i) = (int)(*(frhstmp+i)); thisfield = mxGetFieldByNumber(prhs[6],0,57); BPC_psi = mxGetPr(thisfield); /* BPC_psi */ thisfield = mxGetFieldByNumber(prhs[6],0,58); BPC_phi1 = mxGetPr(thisfield); /* BPC_phi1 */ thisfield = mxGetFieldByNumber(prhs[6],0,59); BPC_phi2 = mxGetPr(thisfield); /* BPC_phi2 */ /* Sparse matrix as returnvalue */ plhs[0] = mxCreateSparse(ncoords+3,ncoords+bfreep+2,ncoords*ncoords,mxREAL); pr = mxGetPr(plhs[0]); ir = mxGetIr(plhs[0]); jc = mxGetJc(plhs[0]); *jc = 0; /* Parameters for rhs-evaluation-call to Matlab */ evalrhs[0] = (struct mxArray_tag*) prhs[0]; evalrhs[1] = mxCreateDoubleMatrix(1,1,mxREAL); zero = mxGetPr(evalrhs[1]); *zero = 0; /*evalrhs[2] = mxCreateDoubleMatrix(nphase+ActiveParams,1,mxREAL);*/ evalrhs[2] = mxCreateDoubleMatrix(nphase,1,mxREAL); xtmp = mxGetPr(evalrhs[2]); for (i=0; i<mxGetNumberOfElements(prhs[2]); i++) { evalrhs[i+3] = mxCreateDoubleMatrix(1,1,mxREAL); ptmp = mxGetPr(evalrhs[i+3]); *ptmp = *(p+i); } /* Parameters for jacobian-call to Matlab */ jacrhs[0] = (struct mxArray_tag*) prhs[0]; jacrhs[1] = (struct mxArray_tag*) prhs[8]; jacrhs[2] = evalrhs[2]; jacrhs[3] = (struct mxArray_tag*) prhs[7]; jacrhs[4] = mxCreateDoubleMatrix(nfreep,1,mxREAL); xtmpn = mxGetPr(jacrhs[4]); for (i=0; i<nfreep; i++) *(xtmpn+i) = *(ActiveParams+i); /* Parameters for jacobianBR-call to Matlab */ jacrhsBR[0] = (struct mxArray_tag*) prhs[0]; jacrhsBR[1] = (struct mxArray_tag*) prhs[10]; jacrhsBR[2] = evalrhs[2]; jacrhsBR[3] = (struct mxArray_tag*) prhs[7]; jacrhsBR[4] = mxCreateDoubleMatrix(nfreep,1,mxREAL); xtmpn = mxGetPr(jacrhsBR[4]); for (i=0; i<nfreep; i++) *(xtmpn+i) = *(ActiveParams+i); jacrhsBR[5] = mxCreateDoubleMatrix(bfreep,1,mxREAL); xtmpm = mxGetPr(jacrhsBR[5]); for (i=0; i<bfreep; i++) *(xtmpm+i) = *(BranchParam+i); filled = 0; /* Help-variable that will be used in storage-procedure */ elementcounter = 0; /* Counts number of elements already stored in sparse matrix */ tmpperiod = *(mxGetPr(prhs[3])); /* Other memory allocations */ /* ------------------------ */ range1 = mxCalloc(ncol+1,sizeof(int)); range2 = mxCalloc(ncol*nphase,sizeof(int)); range3 = mxCalloc((ncol+1)*nphase,sizeof(int)); range4 = mxCalloc(nphase,sizeof(int)); tempmatrix = mxCalloc(nphase*nphase*ncol,sizeof(double)); jac = mxCalloc(nphase*nphase,sizeof(double)); if (nfreep == 1) { jacp = mxCalloc(nphase*bfreep,sizeof(double)); sysjacp = mxCalloc(nphase*ncol*bfreep,sizeof(double)); freepcols = mxCalloc((tps-1)*bfreep*nphase+1,sizeof(double)); } else { jacp = mxCalloc(nphase*nfreep,sizeof(double)); sysjacp = mxCalloc(nphase*ncol*nfreep,sizeof(double)); freepcols = mxCalloc((tps-1)*nfreep*nphase+1,sizeof(double)); } icjac = mxCalloc(ncoords,sizeof(double)); frhs = mxCalloc(nphase*ncol,sizeof(double)); sysjac = mxCalloc(nphase*ncol*(ncol+1)*nphase,sizeof(double)); Tcol = mxCalloc((tps-1)*nphase+1,sizeof(double)); /* Compute third component: the integral constraint */ /* ------------------------------------------------ */ /* Storage in sparse matrix is done later on */ /* Define some ranges */ for (i=0; i<(ncol+1); i++) { *(range1+i) = i; for (j=0; j<nphase; j++) *(range3+i*nphase+j) = i*nphase+j; } for (i=0; i<ntst; i++) { /* Compute elements of third component */ for (j=0; j<(ncol+1)*nphase; j++) *(icjac + *(range3+j)) = *(icjac + *(range3+j)) + (*(dt+i)) * (*(upoldp+(*range1)*nphase+j)) * (*(pwi+j)); /* Shift the ranges to next intervals */ for (j=0; j<ncol+1; j++) { *(range1+j) = *(range1+j) + ncol; for (k=0; k<nphase; k++) *(range3+j*nphase+k) = *(range3+j*nphase+k) + ncol*nphase; } } /* Compute first component */ /* ----------------------- */ /* Define some ranges*/ for (i=0; i<(ncol+1); i++) { *(range1+i) = i; if (i < ncol) for (j=0; j<nphase; j++) { *(range2+i*nphase+j) = i*nphase+j; *(range3+i*nphase+j) = i*nphase+j; } else for (j=0; j<nphase; j++) *(range3+i*nphase+j) = i*nphase+j; } /* Actual computation of component elements */ for (i=0; i<nt