Author: Chris Douglas (@cmdoug) christopher.douglas@duke.edu
ff-mpirun -np 4 foldcontinue.md -param <PARAM1> -param2 <PARAM2> -fi <FILEIN> -fo <FILEOUT>ff-mpirun -np 4 foldcontinue.md -param <PARAM1> -param2 <PARAM2> -fi <FILEIN> -fo <FILEOUT> -mo <MESHOUT>NOTE: This file should not be changed unless you know what you're doing.
SEE ALSO: modecompute.md, basecontinue.md, foldcompute.md, hopfcontinue.md, fohocompute.md
load "iovtk"
load "PETSc"
include "settings.idp"
include "macros_bifbox.idp"
// arguments
string meshin = getARGV("-mi", ""); // input meshfile with extension
string meshout = getARGV("-mo", "");
string filein = getARGV("-fi", "");
string fileout = getARGV("-fo", filein);
int contorder = getARGV("-contorder", 1);
bool normalform = getARGV("-nf", 1);
int count = getARGV("-count", 0);
int savecount = getARGV("-scount", 1);
int maxcount = getARGV("-maxcount", 100);
real h0 = getARGV("-h0", 1.0);
string param = getARGV("-param", "");
string param2 = getARGV("-param2", "");
string adaptto = getARGV("-adaptto","b");
real fmax = getARGV("-fmax", 2.0);
real kappamax = getARGV("-kmax", 1.0);
real deltamax = getARGV("-dmax", 4.0);
real anglemax = getARGV("-amax", 30.)*pi/180.0;
real monotone = getARGV("-mono", 0.0);
real eps = getARGV("-eps", 1e-7);
real eps2 = getARGV("-eps2", 1e-7);
bool stricttangent = bool(getARGV("-stricttangent", 1));
int snesmaxit = getARGV("-snes_max_it", 10);
string sneslinesearchtype = getARGV("-snes_linesearch_type","basic");
real[string] alpha;
real beta;
real paramtarget = getARGV("-paramtarget",1.0);
real param2target = getARGV("-param2target",1.0);
real TGV = getARGV("-tgv", -1);
bool stopflag = false;
bool forcesave = false;
// Load mesh, make FE basis
string fileroot, fileext = parsefilename(filein, fileroot); //extract file name and extension
parsefilename(fileout, fileout); // trim extension from output file, if given
if(meshin == "") meshin = readmeshname(workdir + filein); // get mesh file
string meshroot, meshext = parsefilename(meshin, meshroot);
parsefilename(meshout, meshroot); // trim extension from output mesh, if given
if(count > 0) {
fileroot = fileroot(0:fileroot.rfind("_" + count)-1); // get file root
meshroot = meshroot(0:meshroot.rfind("_" + count)-1); // get file root
}
assert(fileext == "fold" || fileext == "foho");
Th = readmeshN(workdir + meshin);
Thg = Th;
DmeshCreate(Th);
restu = restrict(XMh, XMhg, n2o);
XMh defu(ub), defu(um), defu(uma), defu(um2), defu(um3);
if (count == 0){
if(fileext == "fold"){
ub[] = loadfold(fileroot, meshin, um[], uma[], alpha, beta);
}
else if(fileext == "foho") {
real omega, gamma22, gamma23, beta23;
complex[string] alpha1;
complex beta1, gamma12, gamma13;
complex[int] q1m, q1ma;
ub[] = loadfoho(fileroot, meshin, q1m, q1ma, um[], uma[], sym, omega, alpha1, alpha, beta1, beta, beta23, gamma12, gamma13, gamma22, gamma23);
}
savefold(filein, (savecount > 0 ? fileout : ""), meshin, alpha, beta, false, false);
}
else {
ub[] = loadfold(fileroot + "_" + count, meshin, um[], uma[], alpha, beta);
}
real[int] paramvals(2);
paramvals(0) = getparam(param);
paramvals(1) = getparam(param2);
real paramdiff1 = paramvals(0) - paramtarget;
real paramdiff2 = paramvals(1) - param2target;
// Create distributed Mat
Mat J;
createMatu(Th, J, Pk);
sym = 0;
real[int] ik(sym.n), ik2(sym.n), ik3(sym.n);
real iomega = 0.0, iomega2 = 0.0, iomega3 = 0.0;
include "eqns.idp"
bool adaptflag, adapt = false;
if(meshout != "") adapt = true; // if output meshfile is given, adapt mesh
meshout = meshin; // if no adaptation
// Build bordered block matrix from only Mat components
Mat JlPM(J.n, mpirank == 0 ? 1 : 0), gqPM(J.n, mpirank == 0 ? 1 : 0), glPM(mpirank == 0 ? 1 : 0, mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
Mat JlPMa(J.n + (mpirank == 0 ? 1 : 0), mpirank == 0 ? 1 : 0), yqPMa(J.n + (mpirank == 0 ? 1 : 0), mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
Mat Ja = [[J, JlPM], [gqPM', glPM]], Jaa = [[Ja, JlPMa], [yqPMa', -1.0]]; // make dummy Jacobian
real[int] R(ub[].n), qm(J.n), qma(J.n), pP(J.n), qP(J.n), yqP(Ja.n), yqP0(Ja.n);
int ret, it = 0;
real f, kappa, cosalpha, res, delta, maxdelta, alpha0, beta0;
// FUNCTIONS
func PetscScalar[int] funcRa(PetscScalar[int]& qa) {
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true); // PETSc to FreeFEM
if(mpirank == 0) paramvals = qa(J.n:Jaa.n-1);
broadcast(processor(0), paramvals);
updateparam(param, paramvals(0));
updateparam(param2, paramvals(1));
R = vR(0, XMh, tgv = TGV);
PetscScalar[int] Ra;
ChangeNumbering(J, R, Ra); // FreeFEM to PETSc
J = vJ(XMh, XMh, tgv = -2);
KSPSolve(J, pP, qm);
KSPSolveTranspose(J, qP, qma);
PetscScalar ginv, ginvl = (qP'*qm);
mpiAllReduce(ginvl, ginv, mpiCommWorld, mpiSUM);
qm /= ginv;
qma /= ginv;
Ra.resize(Jaa.n); // Append 0 to residual vector on proc 0
if(mpirank == 0) Ra(J.n:Jaa.n-1) = [1.0/ginv, 0.0];
return Ra;
}
func int funcJa(PetscScalar[int]& qa) {
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true); // PETSc to FreeFEM
if(mpirank == 0) paramvals = qa(J.n:Jaa.n-1);
broadcast(processor(0), paramvals);
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
ChangeNumbering(J, uma[], qma, inverse = true);
updateparam(param, paramvals(0) + eps);
updateparam(param2, paramvals(1));
um2[] = vR(0, XMh, tgv = TGV);
um2[] -= R;
um2[] /= eps;
ChangeNumbering(J, um2[], qm); // FreeFEM to PETSc
matrix<PetscScalar> tempPms = [[qm]]; // dense array to sparse matrix
ChangeOperator(JlPM, tempPms, parent = Ja); // send to Mat
um2[] = vJ(0, XMh, tgv = -10);
updateparam(param, paramvals(0));
um3[] = vJ(0, XMh, tgv = -10);
um2[] -= um3[];
tempPms = [[J(uma[], um2[])/eps]];
ChangeOperator(glPM, tempPms, parent = Ja); // send to Mat
J = vH(XMh, XMh, tgv = 0);
MatMultTranspose(J, qma, qm);
tempPms = [[qm]];
ChangeOperator(gqPM, tempPms, parent = Ja); // send to Mat
J = vJ(XMh, XMh, tgv = TGV);
if (contorder > 0) {
updateparam(param2, paramvals(1) + eps2);
um2[] = vR(0, XMh, tgv = TGV);
um2[] -= R;
um2[] /= eps2;
ChangeNumbering(J, um2[], qm); // FreeFEM to PETSc
yqP(0:J.n-1) = qm;
um2[] = vJ(0, XMh, tgv = -10);
updateparam(param2, paramvals(1));
um2[] -= um3[];
PetscScalar gl = J(uma[], um2[])/eps2;
if (mpirank == 0) yqP(Ja.n-1) = gl;
tempPms = [[yqP]]; // dense array to sparse matrix
ChangeOperator(JlPMa, tempPms, parent = Jaa); // send to Mat
KSPSolve(Ja, yqP, yqP);
tempPms = [[yqP]]; // dense array to sparse matrix
ChangeOperator(yqPMa, tempPms, parent = Jaa); // send to Mat
}
return 0;
}
ConvergenceCheck(param + " = " + paramvals(0) + ",\t" + param2 + " = " + paramvals(1));
// set up Mat parameters
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(Jaa, sparams = "-ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition full"
+ " -prefix_push fieldsplit_1_ -ksp_type preonly -pc_type redundant -redundant_pc_type lu -prefix_pop"
+ " -prefix_push fieldsplit_0_ -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition full -prefix_pop"
+ " -prefix_push fieldsplit_0_fieldsplit_1_ -ksp_type preonly -pc_type redundant -redundant_pc_type lu -prefix_pop"
+ " -prefix_push fieldsplit_0_fieldsplit_0_ " + KSPparams + " -prefix_pop", setup = 1);
set(Ja, prefix = "fieldsplit_0_", setup = 1, parent = Jaa);
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO prefix = "fieldsplit_0_fieldsplit_0_", parent = Ja);
// PREDICTOR
real[int] qa(Jaa.n), qa0;
ChangeNumbering(J, ub[], qa0);
ChangeNumbering(J, um[], qm);
ChangeNumbering(J, uma[], qma);
qa0.resize(Jaa.n);
if(mpirank == 0) qa0(J.n:Jaa.n-1) = paramvals;
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
um2[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um2[], qP);
ChangeNumbering(J, um[], qma, inverse = true, exchange = true);
um2[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um2[], pP);
if (contorder > 0) {
R = vR(0, XMh, tgv = TGV);
funcJa(qa0);
}
else {
matrix<PetscScalar> tempPms = [[yqP]];
ChangeOperator(JlPMa, tempPms, parent = Jaa); // send to Mat
ChangeOperator(yqPMa, tempPms, parent = Jaa); // send to Mat
}
yqP0 = yqP;
alpha0 = alpha[paramnames[0]];
beta0 = beta;
while (!stopflag){
qa = qa0;
if (contorder == 0 && mpirank == 0) qa(Jaa.n-1) += h0;
else if (contorder > 0) {
real h, hl = (yqP0'*yqP0);
mpiAllReduce(hl, h, mpiCommWorld, mpiSUM);
h = h0/sqrt(h + 1.0);
qa(0:Ja.n-1) -= (h*yqP0);
if (mpirank == 0) qa(Jaa.n-1) += h;
}
// CORRECTOR LOOP
adaptflag = false;
SNESSolve(Jaa, funcJa, funcRa, qa, convergence = funcConvergence, reason = ret,
sparams = "-snes_linesearch_type " + sneslinesearchtype + " -options_left no -snes_converged_reason -snes_max_it " + snesmaxit); // solve nonlinear problem with SNES
if (ret > 0) {
++count;
if (maxcount > 0) stopflag = (count >= maxcount);
else if ((paramvals(0) - paramtarget)*paramdiff1 <= 0 || (paramvals(1) - param2target)*paramdiff2 <= 0) stopflag = true;
h0 /= f;
if (cosalpha < 0 && contorder > 0) {
h0 *= -1.0;
if(mpirank == 0) cout << "\tOrientation reversed." << endl;
forcesave = true;
}
if (adapt && (count % savecount == 0)){
meshout = meshroot + "_" + count + "." + meshext;
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true);
if(mpirank == 0) paramvals = qa(J.n:Jaa.n-1);
broadcast(processor(0), paramvals);
updateparam(param, paramvals(0));
updateparam(param2, paramvals(1));
ChangeNumbering(J, um[], qm, inverse = true);
ChangeNumbering(J, uma[], qma, inverse = true);
ChangeNumbering(J, um2[], yqP(0:J.n-1), inverse = true);
real yparamval;
if(mpirank == 0) yparamval = yqP(Ja.n-1);
XMhg defu(uG), defu(umG), defu(umaG), defu(yG), defu(tempu), defu(uoG);
tempu[](restu) = ub[]; // populate local portion of global soln
mpiAllReduce(tempu[], uG[], mpiCommWorld, mpiSUM);
tempu[](restu) = um2[]; // populate local portion of global soln
mpiAllReduce(tempu[], yG[], mpiCommWorld, mpiSUM);
tempu[](restu) = um[]; // populate local portion of global soln
mpiAllReduce(tempu[], umG[], mpiCommWorld, mpiSUM); //aggregate local solns into global soln
tempu[](restu) = uma[]; // populate local portion of global soln
mpiAllReduce(tempu[], umaG[], mpiCommWorld, mpiSUM); //aggregate local solns into global soln
if(mpirank == 0) { // Perform mesh adaptation (serially) on processor 0
if(adaptto == "bo") {
defu(tempu) = initu(defu(umaG)'*defu(umaG));
tempu[] = sqrt(tempu[]);
uoG[] = (umG[].*umG[]);
uoG[] = sqrt(uoG[]);
uoG[] .*= tempu[];
}
IFMACRO(dimension,2)
if(adaptto == "b") Thg = adaptmesh(Thg, adaptu(uG), adaptmeshoptions);
else if(adaptto == "by") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), adaptmeshoptions);
else if(adaptto == "bd") Thg = adaptmesh(Thg, adaptu(uG), adaptu(umG), adaptmeshoptions);
else if(adaptto == "ba") Thg = adaptmesh(Thg, adaptu(uG), adaptu(umaG), adaptmeshoptions);
else if(adaptto == "byd") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), adaptu(umG), adaptmeshoptions);
else if(adaptto == "bya") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), adaptu(umaG), adaptmeshoptions);
else if(adaptto == "bda") Thg = adaptmesh(Thg, adaptu(uG), adaptu(umG), adaptu(umaG), adaptmeshoptions);
else if(adaptto == "byda") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), adaptu(umG), adaptu(umaG), adaptmeshoptions);
else if(adaptto == "bo") Thg = adaptmesh(Thg, adaptu(uG), adaptu(uoG), adaptmeshoptions);
ENDIFMACRO
IFMACRO(dimension,3)
//NOTE: 3D mesh adaptation is still under development.
load "mshmet"
load "mmg"
real anisomax = getARGV("-anisomax",1.0);
real[int] met((bool(anisomax > 1) ? 6 : 1)*Thg.nv);
if(adaptto == "b") met = mshmet(Thg, adaptu(uG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
else if(adaptto == "by") met = mshmet(Thg, adaptu(uG), adaptu(yG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
else if(adaptto == "bd") met = mshmet(Thg, adaptu(uG), adaptu(umG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
else if(adaptto == "ba") met = mshmet(Thg, adaptu(uG), adaptu(umaG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
else if(adaptto == "byd") met = mshmet(Thg, adaptu(uG), adaptu(yG), adaptu(umG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
else if(adaptto == "bya") met = mshmet(Thg, adaptu(uG), adaptu(yG), adaptu(umaG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
else if(adaptto == "bda") met = mshmet(Thg, adaptu(uG), adaptu(umG), adaptu(umaG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
else if(adaptto == "byda") met = mshmet(Thg, adaptu(uG), adaptu(yG), adaptu(umG), adaptu(umaG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
else if(adaptto == "bo") met = mshmet(Thg, adaptu(uG), adaptu(uoG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0),hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
if(anisomax > 1.0) {
load "aniso"
boundaniso(6, met, anisomax);
}
Thg = mmg3d(Thg, metric = met, hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), hgrad = -1, verbose = verbosity-(verbosity==0));
ENDIFMACRO
}
broadcast(processor(0), Thg); // broadcast global mesh to all processors
defu(uG) = defu(uG); //interpolate global solution from old mesh to new mesh
defu(yG) = defu(yG); //interpolate global solution from old mesh to new mesh
defu(umG) = defu(umG); //interpolate global solution from old mesh to new mesh
defu(umaG) = defu(umaG); //interpolate global solution from old mesh to new mesh
Th = Thg;
Mat Adapt;
createMatu(Th, Adapt, Pk);
J = Adapt;
defu(ub) = initu(0.0);
defu(um) = initu(0.0);
defu(uma) = initu(0.0);
defu(um2) = initu(0.0);
defu(um3) = initu(0.0);
restu.resize(ub[].n); // Change size of restriction operator
restu = restrict(XMh, XMhg, n2o); // Compute new restriction from global mesh to local mesh
ub[] = uG[](restu);
um2[] = yG[](restu);
um[] = umG[](restu);
uma[] = umaG[](restu);
Mat Adapt1(J.n, mpirank == 0 ? 1 : 0), Adapt2(J.n, mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
Mat Adapt3(J.n + (mpirank == 0 ? 1 : 0), mpirank == 0 ? 1 : 0), Adapt4(J.n + (mpirank == 0 ? 1 : 0), mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
JlPM = Adapt1;
gqPM = Adapt2;
JlPMa = Adapt3;
yqPMa = Adapt4;
Ja = [[J, JlPM], [gqPM', glPM]]; // make dummy Jacobian
Jaa = [[Ja, JlPMa], [yqPMa', -1.0]]; // make dummy Jacobian
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO prefix = "fieldsplit_0_fieldsplit_0_", parent = Ja);
set(Ja, prefix = "fieldsplit_0_", parent = Jaa);
ChangeNumbering(J, ub[], qa);
qa.resize(Jaa.n);
if(mpirank == 0) qa(J.n:Jaa.n-1) = paramvals;
ChangeNumbering(J, um[], qm);
ChangeNumbering(J, uma[], qma);
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true);
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
um3[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um3[], qP);
ChangeNumbering(J, um[], qma, inverse = true, exchange = true);
um3[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um3[], pP);
R.resize(ub[].n);
ChangeNumbering(J, um2[], yqP);
yqP.resize(Ja.n);
if (mpirank==0) yqP(Ja.n-1) = yparamval;
yqP0.resize(Ja.n);
yqP0 = yqP;
qa0.resize(Jaa.n);
if (contorder == 0) {
matrix<PetscScalar> tempPms = [[yqP]];
ChangeOperator(JlPMa, tempPms, parent = Jaa); // send to Mat
ChangeOperator(yqPMa, tempPms, parent = Jaa); // send to Mat
}
adaptflag = true;
SNESSolve(Jaa, funcJa, funcRa, qa, convergence = funcConvergence, reason = ret,
sparams = "-snes_linesearch_type " + sneslinesearchtype + " -options_left no -snes_converged_reason"); // solve nonlinear problem with SNES
assert(ret > 0);
if(mpirank==0) { // Save adapted mesh
cout << " Saving adapted mesh '" + meshout + "' in '" + workdir + "'." << endl;
savemesh(Thg, workdir + meshout);
}
}
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true);
if(mpirank == 0) paramvals = qa(J.n:Jaa.n-1);
broadcast(processor(0), paramvals);
updateparam(param, paramvals(0));
updateparam(param2, paramvals(1));
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
um2[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um[], qm, inverse = true);
ChangeNumbering(J, uma[], qma, inverse = true);
real Mnorm = sqrt(J(um[], um2[]));
um[] /= Mnorm; // so that <um[],M*um[]> = 1
uma[] *= (Mnorm/J(um2[], uma[])); // so that <uma[],M*um[]> = 1
ChangeNumbering(J, um[], qm);
ChangeNumbering(J, uma[], qma);
if(normalform){
// 2nd-order
// A: base modification due to parameter changes
if(paramnames[0] != ""){
for (int k = 0; k < paramnames.n; ++k){
real paramval = getparam(paramnames[k]);
updateparam(paramnames[k], paramval + eps);
um2[] = vR(0, XMh, tgv = TGV);
updateparam(paramnames[k], paramval);
um2[] -= R;
alpha[paramnames[k]] = J(uma[], um2[])/eps;
}
}
// B: base modification due to quadratic nonlinear interaction
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
um2[] = um[];
um3[] = vH(0, XMh, tgv = -10);
beta = 0.5*J(uma[], um3[]);
}
else {
for (int k = 0; k < paramnames.n; ++k){
alpha[paramnames[k]] = 0.0;
}
beta = 0.0;
}
if (beta*beta0 < 0) {
if(mpirank == 0) cout << "\tCusp bifurcation detected." << endl;
forcesave = true;
}
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true);
ChangeNumbering(J, um[], qm, inverse = true);
savefold(fileout + "_" + count + (forcesave ? "specialpt" : ""), fileout, meshout, alpha, beta, ((count % savecount == 0) || forcesave || stopflag), true);
if (stopflag) break;
forcesave = false;
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true);
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
um2[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um2[], qP);
ChangeNumbering(J, um[], qma, inverse = true, exchange = true);
um2[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um2[], pP);
it = 0;
if (stricttangent && contorder > 0) funcJa(qa);
yqP0 = yqP;
qa0 = qa;
alpha0 = alpha[paramnames[0]];
beta0 = beta;
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
}
else h0 /= fmax;
}