1 #include <../src/tao/leastsquares/impls/pounders/pounders.h> 2 3 static PetscErrorCode pounders_h(Tao subtao, Vec v, Mat H, Mat Hpre, void *ctx) 4 { 5 PetscFunctionBegin; 6 PetscFunctionReturn(0); 7 } 8 9 static PetscErrorCode pounders_fg(Tao subtao, Vec x, PetscReal *f, Vec g, void *ctx) 10 { 11 TAO_POUNDERS *mfqP = (TAO_POUNDERS*)ctx; 12 PetscReal d1,d2; 13 14 PetscFunctionBegin; 15 /* g = A*x (add b later)*/ 16 CHKERRQ(MatMult(mfqP->subH,x,g)); 17 18 /* f = 1/2 * x'*(Ax) + b'*x */ 19 CHKERRQ(VecDot(x,g,&d1)); 20 CHKERRQ(VecDot(mfqP->subb,x,&d2)); 21 *f = 0.5 *d1 + d2; 22 23 /* now g = g + b */ 24 CHKERRQ(VecAXPY(g, 1.0, mfqP->subb)); 25 PetscFunctionReturn(0); 26 } 27 28 static PetscErrorCode pounders_feval(Tao tao, Vec x, Vec F, PetscReal *fsum) 29 { 30 TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data; 31 PetscInt i,row,col; 32 PetscReal fr,fc; 33 34 PetscFunctionBegin; 35 CHKERRQ(TaoComputeResidual(tao,x,F)); 36 if (tao->res_weights_v) { 37 CHKERRQ(VecPointwiseMult(mfqP->workfvec,tao->res_weights_v,F)); 38 CHKERRQ(VecDot(mfqP->workfvec,mfqP->workfvec,fsum)); 39 } else if (tao->res_weights_w) { 40 *fsum=0; 41 for (i=0;i<tao->res_weights_n;i++) { 42 row=tao->res_weights_rows[i]; 43 col=tao->res_weights_cols[i]; 44 CHKERRQ(VecGetValues(F,1,&row,&fr)); 45 CHKERRQ(VecGetValues(F,1,&col,&fc)); 46 *fsum += tao->res_weights_w[i]*fc*fr; 47 } 48 } else { 49 CHKERRQ(VecDot(F,F,fsum)); 50 } 51 CHKERRQ(PetscInfo(tao,"Least-squares residual norm: %20.19e\n",(double)*fsum)); 52 PetscCheck(!PetscIsInfOrNanReal(*fsum),PETSC_COMM_SELF,PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode gqtwrap(Tao tao,PetscReal *gnorm, PetscReal *qmin) 57 { 58 #if defined(PETSC_USE_REAL_SINGLE) 59 PetscReal atol=1.0e-5; 60 #else 61 PetscReal atol=1.0e-10; 62 #endif 63 PetscInt info,its; 64 TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data; 65 66 PetscFunctionBegin; 67 if (!mfqP->usegqt) { 68 PetscReal maxval; 69 PetscInt i,j; 70 71 CHKERRQ(VecSetValues(mfqP->subb,mfqP->n,mfqP->indices,mfqP->Gres,INSERT_VALUES)); 72 CHKERRQ(VecAssemblyBegin(mfqP->subb)); 73 CHKERRQ(VecAssemblyEnd(mfqP->subb)); 74 75 CHKERRQ(VecSet(mfqP->subx,0.0)); 76 77 CHKERRQ(VecSet(mfqP->subndel,-1.0)); 78 CHKERRQ(VecSet(mfqP->subpdel,+1.0)); 79 80 /* Complete the lower triangle of the Hessian matrix */ 81 for (i=0;i<mfqP->n;i++) { 82 for (j=i+1;j<mfqP->n;j++) { 83 mfqP->Hres[j+mfqP->n*i] = mfqP->Hres[mfqP->n*j+i]; 84 } 85 } 86 CHKERRQ(MatSetValues(mfqP->subH,mfqP->n,mfqP->indices,mfqP->n,mfqP->indices,mfqP->Hres,INSERT_VALUES)); 87 CHKERRQ(MatAssemblyBegin(mfqP->subH,MAT_FINAL_ASSEMBLY)); 88 CHKERRQ(MatAssemblyEnd(mfqP->subH,MAT_FINAL_ASSEMBLY)); 89 90 CHKERRQ(TaoResetStatistics(mfqP->subtao)); 91 /* CHKERRQ(TaoSetTolerances(mfqP->subtao,*gnorm,*gnorm,PETSC_DEFAULT)); */ 92 /* enforce bound constraints -- experimental */ 93 if (tao->XU && tao->XL) { 94 CHKERRQ(VecCopy(tao->XU,mfqP->subxu)); 95 CHKERRQ(VecAXPY(mfqP->subxu,-1.0,tao->solution)); 96 CHKERRQ(VecScale(mfqP->subxu,1.0/mfqP->delta)); 97 CHKERRQ(VecCopy(tao->XL,mfqP->subxl)); 98 CHKERRQ(VecAXPY(mfqP->subxl,-1.0,tao->solution)); 99 CHKERRQ(VecScale(mfqP->subxl,1.0/mfqP->delta)); 100 101 CHKERRQ(VecPointwiseMin(mfqP->subxu,mfqP->subxu,mfqP->subpdel)); 102 CHKERRQ(VecPointwiseMax(mfqP->subxl,mfqP->subxl,mfqP->subndel)); 103 } else { 104 CHKERRQ(VecCopy(mfqP->subpdel,mfqP->subxu)); 105 CHKERRQ(VecCopy(mfqP->subndel,mfqP->subxl)); 106 } 107 /* Make sure xu > xl */ 108 CHKERRQ(VecCopy(mfqP->subxl,mfqP->subpdel)); 109 CHKERRQ(VecAXPY(mfqP->subpdel,-1.0,mfqP->subxu)); 110 CHKERRQ(VecMax(mfqP->subpdel,NULL,&maxval)); 111 PetscCheck(maxval <= 1e-10,PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_OUTOFRANGE,"upper bound < lower bound in subproblem"); 112 /* Make sure xu > tao->solution > xl */ 113 CHKERRQ(VecCopy(mfqP->subxl,mfqP->subpdel)); 114 CHKERRQ(VecAXPY(mfqP->subpdel,-1.0,mfqP->subx)); 115 CHKERRQ(VecMax(mfqP->subpdel,NULL,&maxval)); 116 PetscCheck(maxval <= 1e-10,PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_OUTOFRANGE,"initial guess < lower bound in subproblem"); 117 118 CHKERRQ(VecCopy(mfqP->subx,mfqP->subpdel)); 119 CHKERRQ(VecAXPY(mfqP->subpdel,-1.0,mfqP->subxu)); 120 CHKERRQ(VecMax(mfqP->subpdel,NULL,&maxval)); 121 PetscCheck(maxval <= 1e-10,PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_OUTOFRANGE,"initial guess > upper bound in subproblem"); 122 123 CHKERRQ(TaoSolve(mfqP->subtao)); 124 CHKERRQ(TaoGetSolutionStatus(mfqP->subtao,NULL,qmin,NULL,NULL,NULL,NULL)); 125 126 /* test bounds post-solution*/ 127 CHKERRQ(VecCopy(mfqP->subxl,mfqP->subpdel)); 128 CHKERRQ(VecAXPY(mfqP->subpdel,-1.0,mfqP->subx)); 129 CHKERRQ(VecMax(mfqP->subpdel,NULL,&maxval)); 130 if (maxval > 1e-5) { 131 CHKERRQ(PetscInfo(tao,"subproblem solution < lower bound\n")); 132 tao->reason = TAO_DIVERGED_TR_REDUCTION; 133 } 134 135 CHKERRQ(VecCopy(mfqP->subx,mfqP->subpdel)); 136 CHKERRQ(VecAXPY(mfqP->subpdel,-1.0,mfqP->subxu)); 137 CHKERRQ(VecMax(mfqP->subpdel,NULL,&maxval)); 138 if (maxval > 1e-5) { 139 CHKERRQ(PetscInfo(tao,"subproblem solution > upper bound\n")); 140 tao->reason = TAO_DIVERGED_TR_REDUCTION; 141 } 142 } else { 143 gqt(mfqP->n,mfqP->Hres,mfqP->n,mfqP->Gres,1.0,mfqP->gqt_rtol,atol,mfqP->gqt_maxits,gnorm,qmin,mfqP->Xsubproblem,&info,&its,mfqP->work,mfqP->work2, mfqP->work3); 144 } 145 *qmin *= -1; 146 PetscFunctionReturn(0); 147 } 148 149 static PetscErrorCode pounders_update_res(Tao tao) 150 { 151 TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data; 152 PetscInt i,row,col; 153 PetscBLASInt blasn=mfqP->n,blasn2=blasn*blasn,blasm=mfqP->m,ione=1; 154 PetscReal zero=0.0,one=1.0,wii,factor; 155 156 PetscFunctionBegin; 157 for (i=0;i<mfqP->n;i++) { 158 mfqP->Gres[i]=0; 159 } 160 for (i=0;i<mfqP->n*mfqP->n;i++) { 161 mfqP->Hres[i]=0; 162 } 163 164 /* Compute Gres= sum_ij[wij * (cjgi + cigj)] */ 165 if (tao->res_weights_v) { 166 /* Vector(diagonal) weights: gres = sum_i(wii*ci*gi) */ 167 for (i=0;i<mfqP->m;i++) { 168 CHKERRQ(VecGetValues(tao->res_weights_v,1,&i,&factor)); 169 factor=factor*mfqP->C[i]; 170 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn,&factor,&mfqP->Fdiff[blasn*i],&ione,mfqP->Gres,&ione)); 171 } 172 173 /* compute Hres = sum_ij [wij * (*ci*Hj + cj*Hi + gi gj' + gj gi') ] */ 174 /* vector(diagonal weights) Hres = sum_i(wii*(ci*Hi + gi * gi')*/ 175 for (i=0;i<mfqP->m;i++) { 176 CHKERRQ(VecGetValues(tao->res_weights_v,1,&i,&wii)); 177 if (tao->niter>1) { 178 factor=wii*mfqP->C[i]; 179 /* add wii * ci * Hi */ 180 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn2,&factor,&mfqP->H[i],&blasm,mfqP->Hres,&ione)); 181 } 182 /* add wii * gi * gi' */ 183 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&blasn,&blasn,&ione,&wii,&mfqP->Fdiff[blasn*i],&blasn,&mfqP->Fdiff[blasn*i],&blasn,&one,mfqP->Hres,&blasn)); 184 } 185 } else if (tao->res_weights_w) { 186 /* General case: .5 * Gres= sum_ij[wij * (cjgi + cigj)] */ 187 for (i=0;i<tao->res_weights_n;i++) { 188 row=tao->res_weights_rows[i]; 189 col=tao->res_weights_cols[i]; 190 191 factor = tao->res_weights_w[i]*mfqP->C[col]/2.0; 192 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn,&factor,&mfqP->Fdiff[blasn*row],&ione,mfqP->Gres,&ione)); 193 factor = tao->res_weights_w[i]*mfqP->C[row]/2.0; 194 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn,&factor,&mfqP->Fdiff[blasn*col],&ione,mfqP->Gres,&ione)); 195 } 196 197 /* compute Hres = sum_ij [wij * (*ci*Hj + cj*Hi + gi gj' + gj gi') ] */ 198 /* .5 * sum_ij [wij * (*ci*Hj + cj*Hi + gi gj' + gj gi') ] */ 199 for (i=0;i<tao->res_weights_n;i++) { 200 row=tao->res_weights_rows[i]; 201 col=tao->res_weights_cols[i]; 202 factor=tao->res_weights_w[i]/2.0; 203 /* add wij * gi gj' + wij * gj gi' */ 204 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&blasn,&blasn,&ione,&factor,&mfqP->Fdiff[blasn*row],&blasn,&mfqP->Fdiff[blasn*col],&blasn,&one,mfqP->Hres,&blasn)); 205 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&blasn,&blasn,&ione,&factor,&mfqP->Fdiff[blasn*col],&blasn,&mfqP->Fdiff[blasn*row],&blasn,&one,mfqP->Hres,&blasn)); 206 } 207 if (tao->niter > 1) { 208 for (i=0;i<tao->res_weights_n;i++) { 209 row=tao->res_weights_rows[i]; 210 col=tao->res_weights_cols[i]; 211 212 /* add wij*cj*Hi */ 213 factor = tao->res_weights_w[i]*mfqP->C[col]/2.0; 214 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn2,&factor,&mfqP->H[row],&blasm,mfqP->Hres,&ione)); 215 216 /* add wij*ci*Hj */ 217 factor = tao->res_weights_w[i]*mfqP->C[row]/2.0; 218 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn2,&factor,&mfqP->H[col],&blasm,mfqP->Hres,&ione)); 219 } 220 } 221 } else { 222 /* Default: Gres= sum_i[cigi] = G*c' */ 223 CHKERRQ(PetscInfo(tao,"Identity weights\n")); 224 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&blasn,&blasm,&one,mfqP->Fdiff,&blasn,mfqP->C,&ione,&zero,mfqP->Gres,&ione)); 225 226 /* compute Hres = sum_ij [wij * (*ci*Hj + cj*Hi + gi gj' + gj gi') ] */ 227 /* Hres = G*G' + 0.5 sum {F(xkin,i)*H(:,:,i)} */ 228 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&blasn,&blasn,&blasm,&one,mfqP->Fdiff, &blasn,mfqP->Fdiff, &blasn,&zero,mfqP->Hres,&blasn)); 229 230 /* sum(F(xkin,i)*H(:,:,i)) */ 231 if (tao->niter>1) { 232 for (i=0;i<mfqP->m;i++) { 233 factor = mfqP->C[i]; 234 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&blasn2,&factor,&mfqP->H[i],&blasm,mfqP->Hres,&ione)); 235 } 236 } 237 } 238 PetscFunctionReturn(0); 239 } 240 241 static PetscErrorCode phi2eval(PetscReal *x, PetscInt n, PetscReal *phi) 242 { 243 /* Phi = .5*[x(1)^2 sqrt(2)*x(1)*x(2) ... sqrt(2)*x(1)*x(n) ... x(2)^2 sqrt(2)*x(2)*x(3) .. x(n)^2] */ 244 PetscInt i,j,k; 245 PetscReal sqrt2 = PetscSqrtReal(2.0); 246 247 PetscFunctionBegin; 248 j=0; 249 for (i=0;i<n;i++) { 250 phi[j] = 0.5 * x[i]*x[i]; 251 j++; 252 for (k=i+1;k<n;k++) { 253 phi[j] = x[i]*x[k]/sqrt2; 254 j++; 255 } 256 } 257 PetscFunctionReturn(0); 258 } 259 260 static PetscErrorCode getquadpounders(TAO_POUNDERS *mfqP) 261 { 262 /* Computes the parameters of the quadratic Q(x) = c + g'*x + 0.5*x*G*x' 263 that satisfies the interpolation conditions Q(X[:,j]) = f(j) 264 for j=1,...,m and with a Hessian matrix of least Frobenius norm */ 265 266 /* NB --we are ignoring c */ 267 PetscInt i,j,k,num,np = mfqP->nmodelpoints; 268 PetscReal one = 1.0,zero=0.0,negone=-1.0; 269 PetscBLASInt blasnpmax = mfqP->npmax; 270 PetscBLASInt blasnplus1 = mfqP->n+1; 271 PetscBLASInt blasnp = np; 272 PetscBLASInt blasint = mfqP->n*(mfqP->n+1) / 2; 273 PetscBLASInt blasint2 = np - mfqP->n-1; 274 PetscBLASInt info,ione=1; 275 PetscReal sqrt2 = PetscSqrtReal(2.0); 276 277 PetscFunctionBegin; 278 for (i=0;i<mfqP->n*mfqP->m;i++) { 279 mfqP->Gdel[i] = 0; 280 } 281 for (i=0;i<mfqP->n*mfqP->n*mfqP->m;i++) { 282 mfqP->Hdel[i] = 0; 283 } 284 285 /* factor M */ 286 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&blasnplus1,&blasnp,mfqP->M,&blasnplus1,mfqP->npmaxiwork,&info)); 287 PetscCheck(info == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACK routine getrf returned with value %d",info); 288 289 if (np == mfqP->n+1) { 290 for (i=0;i<mfqP->npmax-mfqP->n-1;i++) { 291 mfqP->omega[i]=0.0; 292 } 293 for (i=0;i<mfqP->n*(mfqP->n+1)/2;i++) { 294 mfqP->beta[i]=0.0; 295 } 296 } else { 297 /* Let Ltmp = (L'*L) */ 298 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&blasint2,&blasint2,&blasint,&one,&mfqP->L[(mfqP->n+1)*blasint],&blasint,&mfqP->L[(mfqP->n+1)*blasint],&blasint,&zero,mfqP->L_tmp,&blasint)); 299 300 /* factor Ltmp */ 301 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&blasint2,mfqP->L_tmp,&blasint,&info)); 302 PetscCheck(info == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACK routine potrf returned with value %d",info); 303 } 304 305 for (k=0;k<mfqP->m;k++) { 306 if (np != mfqP->n+1) { 307 /* Solve L'*L*Omega = Z' * RESk*/ 308 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&blasnp,&blasint2,&one,mfqP->Z,&blasnpmax,&mfqP->RES[mfqP->npmax*k],&ione,&zero,mfqP->omega,&ione)); 309 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&blasint2,&ione,mfqP->L_tmp,&blasint,mfqP->omega,&blasint2,&info)); 310 PetscCheck(info == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACK routine potrs returned with value %d",info); 311 312 /* Beta = L*Omega */ 313 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&blasint,&blasint2,&one,&mfqP->L[(mfqP->n+1)*blasint],&blasint,mfqP->omega,&ione,&zero,mfqP->beta,&ione)); 314 } 315 316 /* solve M'*Alpha = RESk - N'*Beta */ 317 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&blasint,&blasnp,&negone,mfqP->N,&blasint,mfqP->beta,&ione,&one,&mfqP->RES[mfqP->npmax*k],&ione)); 318 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&blasnplus1,&ione,mfqP->M,&blasnplus1,mfqP->npmaxiwork,&mfqP->RES[mfqP->npmax*k],&blasnplus1,&info)); 319 PetscCheck(info == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACK routine getrs returned with value %d",info); 320 321 /* Gdel(:,k) = Alpha(2:n+1) */ 322 for (i=0;i<mfqP->n;i++) { 323 mfqP->Gdel[i + mfqP->n*k] = mfqP->RES[mfqP->npmax*k + i+1]; 324 } 325 326 /* Set Hdels */ 327 num=0; 328 for (i=0;i<mfqP->n;i++) { 329 /* H[i,i,k] = Beta(num) */ 330 mfqP->Hdel[(i*mfqP->n + i)*mfqP->m + k] = mfqP->beta[num]; 331 num++; 332 for (j=i+1;j<mfqP->n;j++) { 333 /* H[i,j,k] = H[j,i,k] = Beta(num)/sqrt(2) */ 334 mfqP->Hdel[(j*mfqP->n + i)*mfqP->m + k] = mfqP->beta[num]/sqrt2; 335 mfqP->Hdel[(i*mfqP->n + j)*mfqP->m + k] = mfqP->beta[num]/sqrt2; 336 num++; 337 } 338 } 339 } 340 PetscFunctionReturn(0); 341 } 342 343 static PetscErrorCode morepoints(TAO_POUNDERS *mfqP) 344 { 345 /* Assumes mfqP->model_indices[0] is minimum index 346 Finishes adding points to mfqP->model_indices (up to npmax) 347 Computes L,Z,M,N 348 np is actual number of points in model (should equal npmax?) */ 349 PetscInt point,i,j,offset; 350 PetscInt reject; 351 PetscBLASInt blasn=mfqP->n,blasnpmax=mfqP->npmax,blasnplus1=mfqP->n+1,info,blasnmax=mfqP->nmax,blasint,blasint2,blasnp,blasmaxmn; 352 const PetscReal *x; 353 PetscReal normd; 354 355 PetscFunctionBegin; 356 /* Initialize M,N */ 357 for (i=0;i<mfqP->n+1;i++) { 358 CHKERRQ(VecGetArrayRead(mfqP->Xhist[mfqP->model_indices[i]],&x)); 359 mfqP->M[(mfqP->n+1)*i] = 1.0; 360 for (j=0;j<mfqP->n;j++) { 361 mfqP->M[j+1+((mfqP->n+1)*i)] = (x[j] - mfqP->xmin[j]) / mfqP->delta; 362 } 363 CHKERRQ(VecRestoreArrayRead(mfqP->Xhist[mfqP->model_indices[i]],&x)); 364 CHKERRQ(phi2eval(&mfqP->M[1+((mfqP->n+1)*i)],mfqP->n,&mfqP->N[mfqP->n*(mfqP->n+1)/2 * i])); 365 } 366 367 /* Now we add points until we have npmax starting with the most recent ones */ 368 point = mfqP->nHist-1; 369 mfqP->nmodelpoints = mfqP->n+1; 370 while (mfqP->nmodelpoints < mfqP->npmax && point>=0) { 371 /* Reject any points already in the model */ 372 reject = 0; 373 for (j=0;j<mfqP->n+1;j++) { 374 if (point == mfqP->model_indices[j]) { 375 reject = 1; 376 break; 377 } 378 } 379 380 /* Reject if norm(d) >c2 */ 381 if (!reject) { 382 CHKERRQ(VecCopy(mfqP->Xhist[point],mfqP->workxvec)); 383 CHKERRQ(VecAXPY(mfqP->workxvec,-1.0,mfqP->Xhist[mfqP->minindex])); 384 CHKERRQ(VecNorm(mfqP->workxvec,NORM_2,&normd)); 385 normd /= mfqP->delta; 386 if (normd > mfqP->c2) { 387 reject =1; 388 } 389 } 390 if (reject) { 391 point--; 392 continue; 393 } 394 395 CHKERRQ(VecGetArrayRead(mfqP->Xhist[point],&x)); 396 mfqP->M[(mfqP->n+1)*mfqP->nmodelpoints] = 1.0; 397 for (j=0;j<mfqP->n;j++) { 398 mfqP->M[j+1+((mfqP->n+1)*mfqP->nmodelpoints)] = (x[j] - mfqP->xmin[j]) / mfqP->delta; 399 } 400 CHKERRQ(VecRestoreArrayRead(mfqP->Xhist[point],&x)); 401 CHKERRQ(phi2eval(&mfqP->M[1+(mfqP->n+1)*mfqP->nmodelpoints],mfqP->n,&mfqP->N[mfqP->n*(mfqP->n+1)/2 * (mfqP->nmodelpoints)])); 402 403 /* Update QR factorization */ 404 /* Copy M' to Q_tmp */ 405 for (i=0;i<mfqP->n+1;i++) { 406 for (j=0;j<mfqP->npmax;j++) { 407 mfqP->Q_tmp[j+mfqP->npmax*i] = mfqP->M[i+(mfqP->n+1)*j]; 408 } 409 } 410 blasnp = mfqP->nmodelpoints+1; 411 /* Q_tmp,R = qr(M') */ 412 blasmaxmn=PetscMax(mfqP->m,mfqP->n+1); 413 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&blasnp,&blasnplus1,mfqP->Q_tmp,&blasnpmax,mfqP->tau_tmp,mfqP->mwork,&blasmaxmn,&info)); 414 PetscCheck(info == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACK routine geqrf returned with value %d",info); 415 416 /* Reject if min(svd(N*Q(:,n+2:np+1)) <= theta2 */ 417 /* L = N*Qtmp */ 418 blasint2 = mfqP->n * (mfqP->n+1) / 2; 419 /* Copy N to L_tmp */ 420 for (i=0;i<mfqP->n*(mfqP->n+1)/2 * mfqP->npmax;i++) { 421 mfqP->L_tmp[i]= mfqP->N[i]; 422 } 423 /* Copy L_save to L_tmp */ 424 425 /* L_tmp = N*Qtmp' */ 426 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("R","N",&blasint2,&blasnp,&blasnplus1,mfqP->Q_tmp,&blasnpmax,mfqP->tau_tmp,mfqP->L_tmp,&blasint2,mfqP->npmaxwork,&blasnmax,&info)); 427 PetscCheck(info == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACK routine ormqr returned with value %d",info); 428 429 /* Copy L_tmp to L_save */ 430 for (i=0;i<mfqP->npmax * mfqP->n*(mfqP->n+1)/2;i++) { 431 mfqP->L_save[i] = mfqP->L_tmp[i]; 432 } 433 434 /* Get svd for L_tmp(:,n+2:np+1) (L_tmp is modified in process) */ 435 blasint = mfqP->nmodelpoints - mfqP->n; 436 CHKERRQ(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 437 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&blasint2,&blasint,&mfqP->L_tmp[(mfqP->n+1)*blasint2],&blasint2,mfqP->beta,mfqP->work,&blasn,mfqP->work,&blasn,mfqP->npmaxwork,&blasnmax,&info)); 438 CHKERRQ(PetscFPTrapPop()); 439 PetscCheck(info == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACK routine gesvd returned with value %d",info); 440 441 if (mfqP->beta[PetscMin(blasint,blasint2)-1] > mfqP->theta2) { 442 /* accept point */ 443 mfqP->model_indices[mfqP->nmodelpoints] = point; 444 /* Copy Q_tmp to Q */ 445 for (i=0;i<mfqP->npmax* mfqP->npmax;i++) { 446 mfqP->Q[i] = mfqP->Q_tmp[i]; 447 } 448 for (i=0;i<mfqP->npmax;i++) { 449 mfqP->tau[i] = mfqP->tau_tmp[i]; 450 } 451 mfqP->nmodelpoints++; 452 blasnp = mfqP->nmodelpoints; 453 454 /* Copy L_save to L */ 455 for (i=0;i<mfqP->npmax * mfqP->n*(mfqP->n+1)/2;i++) { 456 mfqP->L[i] = mfqP->L_save[i]; 457 } 458 } 459 point--; 460 } 461 462 blasnp = mfqP->nmodelpoints; 463 /* Copy Q(:,n+2:np) to Z */ 464 /* First set Q_tmp to I */ 465 for (i=0;i<mfqP->npmax*mfqP->npmax;i++) { 466 mfqP->Q_tmp[i] = 0.0; 467 } 468 for (i=0;i<mfqP->npmax;i++) { 469 mfqP->Q_tmp[i + mfqP->npmax*i] = 1.0; 470 } 471 472 /* Q_tmp = I * Q */ 473 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("R","N",&blasnp,&blasnp,&blasnplus1,mfqP->Q,&blasnpmax,mfqP->tau,mfqP->Q_tmp,&blasnpmax,mfqP->npmaxwork,&blasnmax,&info)); 474 PetscCheck(info == 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"LAPACK routine ormqr returned with value %d",info); 475 476 /* Copy Q_tmp(:,n+2:np) to Z) */ 477 offset = mfqP->npmax * (mfqP->n+1); 478 for (i=offset;i<mfqP->npmax*mfqP->npmax;i++) { 479 mfqP->Z[i-offset] = mfqP->Q_tmp[i]; 480 } 481 482 if (mfqP->nmodelpoints == mfqP->n + 1) { 483 /* Set L to I_{n+1} */ 484 for (i=0;i<mfqP->npmax * mfqP->n*(mfqP->n+1)/2;i++) { 485 mfqP->L[i] = 0.0; 486 } 487 for (i=0;i<mfqP->n;i++) { 488 mfqP->L[(mfqP->n*(mfqP->n+1)/2)*i + i] = 1.0; 489 } 490 } 491 PetscFunctionReturn(0); 492 } 493 494 /* Only call from modelimprove, addpoint() needs ->Q_tmp and ->work to be set */ 495 static PetscErrorCode addpoint(Tao tao, TAO_POUNDERS *mfqP, PetscInt index) 496 { 497 PetscFunctionBegin; 498 /* Create new vector in history: X[newidx] = X[mfqP->index] + delta*X[index]*/ 499 CHKERRQ(VecDuplicate(mfqP->Xhist[0],&mfqP->Xhist[mfqP->nHist])); 500 CHKERRQ(VecSetValues(mfqP->Xhist[mfqP->nHist],mfqP->n,mfqP->indices,&mfqP->Q_tmp[index*mfqP->npmax],INSERT_VALUES)); 501 CHKERRQ(VecAssemblyBegin(mfqP->Xhist[mfqP->nHist])); 502 CHKERRQ(VecAssemblyEnd(mfqP->Xhist[mfqP->nHist])); 503 CHKERRQ(VecAYPX(mfqP->Xhist[mfqP->nHist],mfqP->delta,mfqP->Xhist[mfqP->minindex])); 504 505 /* Project into feasible region */ 506 if (tao->XU && tao->XL) { 507 CHKERRQ(VecMedian(mfqP->Xhist[mfqP->nHist], tao->XL, tao->XU, mfqP->Xhist[mfqP->nHist])); 508 } 509 510 /* Compute value of new vector */ 511 CHKERRQ(VecDuplicate(mfqP->Fhist[0],&mfqP->Fhist[mfqP->nHist])); 512 CHKMEMQ; 513 CHKERRQ(pounders_feval(tao,mfqP->Xhist[mfqP->nHist],mfqP->Fhist[mfqP->nHist],&mfqP->Fres[mfqP->nHist])); 514 515 /* Add new vector to model */ 516 mfqP->model_indices[mfqP->nmodelpoints] = mfqP->nHist; 517 mfqP->nmodelpoints++; 518 mfqP->nHist++; 519 PetscFunctionReturn(0); 520 } 521 522 static PetscErrorCode modelimprove(Tao tao, TAO_POUNDERS *mfqP, PetscInt addallpoints) 523 { 524 /* modeld = Q(:,np+1:n)' */ 525 PetscInt i,j,minindex=0; 526 PetscReal dp,half=0.5,one=1.0,minvalue=PETSC_INFINITY; 527 PetscBLASInt blasn=mfqP->n, blasnpmax = mfqP->npmax, blask,info; 528 PetscBLASInt blas1=1,blasnmax = mfqP->nmax; 529 530 PetscFunctionBegin; 531 blask = mfqP->nmodelpoints; 532 /* Qtmp = I(n x n) */ 533 for (i=0;i<mfqP->n;i++) { 534 for (j=0;j<mfqP->n;j++) { 535 mfqP->Q_tmp[i + mfqP->npmax*j] = 0.0; 536 } 537 } 538 for (j=0;j<mfqP->n;j++) { 539 mfqP->Q_tmp[j + mfqP->npmax*j] = 1.0; 540 } 541 542 /* Qtmp = Q * I */ 543 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("R","N",&blasn,&blasn,&blask,mfqP->Q,&blasnpmax,mfqP->tau, mfqP->Q_tmp, &blasnpmax, mfqP->npmaxwork,&blasnmax, &info)); 544 545 for (i=mfqP->nmodelpoints;i<mfqP->n;i++) { 546 PetscStackCallBLAS("BLASdot",dp = BLASdot_(&blasn,&mfqP->Q_tmp[i*mfqP->npmax],&blas1,mfqP->Gres,&blas1)); 547 if (dp>0.0) { /* Model says use the other direction! */ 548 for (j=0;j<mfqP->n;j++) { 549 mfqP->Q_tmp[i*mfqP->npmax+j] *= -1; 550 } 551 } 552 /* mfqP->work[i] = Cres+Modeld(i,:)*(Gres+.5*Hres*Modeld(i,:)') */ 553 for (j=0;j<mfqP->n;j++) { 554 mfqP->work2[j] = mfqP->Gres[j]; 555 } 556 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&blasn,&blasn,&half,mfqP->Hres,&blasn,&mfqP->Q_tmp[i*mfqP->npmax], &blas1, &one, mfqP->work2,&blas1)); 557 PetscStackCallBLAS("BLASdot",mfqP->work[i] = BLASdot_(&blasn,&mfqP->Q_tmp[i*mfqP->npmax],&blas1,mfqP->work2,&blas1)); 558 if (i==mfqP->nmodelpoints || mfqP->work[i] < minvalue) { 559 minindex=i; 560 minvalue = mfqP->work[i]; 561 } 562 if (addallpoints != 0) { 563 CHKERRQ(addpoint(tao,mfqP,i)); 564 } 565 } 566 if (!addallpoints) { 567 CHKERRQ(addpoint(tao,mfqP,minindex)); 568 } 569 PetscFunctionReturn(0); 570 } 571 572 static PetscErrorCode affpoints(TAO_POUNDERS *mfqP, PetscReal *xmin,PetscReal c) 573 { 574 PetscInt i,j; 575 PetscBLASInt blasm=mfqP->m,blasj,blask,blasn=mfqP->n,ione=1,info; 576 PetscBLASInt blasnpmax = mfqP->npmax,blasmaxmn; 577 PetscReal proj,normd; 578 const PetscReal *x; 579 580 PetscFunctionBegin; 581 for (i=mfqP->nHist-1;i>=0;i--) { 582 CHKERRQ(VecGetArrayRead(mfqP->Xhist[i],&x)); 583 for (j=0;j<mfqP->n;j++) { 584 mfqP->work[j] = (x[j] - xmin[j])/mfqP->delta; 585 } 586 CHKERRQ(VecRestoreArrayRead(mfqP->Xhist[i],&x)); 587 PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn,mfqP->work,&ione,mfqP->work2,&ione)); 588 PetscStackCallBLAS("BLASnrm2",normd = BLASnrm2_(&blasn,mfqP->work,&ione)); 589 if (normd <= c) { 590 blasj=PetscMax((mfqP->n - mfqP->nmodelpoints),0); 591 if (!mfqP->q_is_I) { 592 /* project D onto null */ 593 blask=(mfqP->nmodelpoints); 594 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("R","N",&ione,&blasn,&blask,mfqP->Q,&blasnpmax,mfqP->tau,mfqP->work2,&ione,mfqP->mwork,&blasm,&info)); 595 PetscCheck(info >= 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"ormqr returned value %d",info); 596 } 597 PetscStackCallBLAS("BLASnrm2",proj = BLASnrm2_(&blasj,&mfqP->work2[mfqP->nmodelpoints],&ione)); 598 599 if (proj >= mfqP->theta1) { /* add this index to model */ 600 mfqP->model_indices[mfqP->nmodelpoints]=i; 601 mfqP->nmodelpoints++; 602 PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasn,mfqP->work,&ione,&mfqP->Q_tmp[mfqP->npmax*(mfqP->nmodelpoints-1)],&ione)); 603 blask=mfqP->npmax*(mfqP->nmodelpoints); 604 PetscStackCallBLAS("BLAScopy",BLAScopy_(&blask,mfqP->Q_tmp,&ione,mfqP->Q,&ione)); 605 blask = mfqP->nmodelpoints; 606 blasmaxmn = PetscMax(mfqP->m,mfqP->n); 607 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&blasn,&blask,mfqP->Q,&blasnpmax,mfqP->tau,mfqP->mwork,&blasmaxmn,&info)); 608 PetscCheck(info >= 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"geqrf returned value %d",info); 609 mfqP->q_is_I = 0; 610 } 611 if (mfqP->nmodelpoints == mfqP->n) { 612 break; 613 } 614 } 615 } 616 617 PetscFunctionReturn(0); 618 } 619 620 static PetscErrorCode TaoSolve_POUNDERS(Tao tao) 621 { 622 TAO_POUNDERS *mfqP = (TAO_POUNDERS *)tao->data; 623 PetscInt i,ii,j,k,l; 624 PetscReal step=1.0; 625 PetscInt low,high; 626 PetscReal minnorm; 627 PetscReal *x,*f; 628 const PetscReal *xmint,*fmin; 629 PetscReal deltaold; 630 PetscReal gnorm; 631 PetscBLASInt info,ione=1,iblas; 632 PetscBool valid,same; 633 PetscReal mdec, rho, normxsp; 634 PetscReal one=1.0,zero=0.0,ratio; 635 PetscBLASInt blasm,blasn,blasncopy,blasnpmax; 636 PetscErrorCode ierr; 637 static PetscBool set = PETSC_FALSE; 638 639 /* n = # of parameters 640 m = dimension (components) of function */ 641 PetscFunctionBegin; 642 ierr = PetscCitationsRegister("@article{UNEDF0,\n" 643 "title = {Nuclear energy density optimization},\n" 644 "author = {Kortelainen, M. and Lesinski, T. and Mor\'e, J. and Nazarewicz, W.\n" 645 " and Sarich, J. and Schunck, N. and Stoitsov, M. V. and Wild, S. },\n" 646 "journal = {Phys. Rev. C},\n" 647 "volume = {82},\n" 648 "number = {2},\n" 649 "pages = {024313},\n" 650 "numpages = {18},\n" 651 "year = {2010},\n" 652 "month = {Aug},\n" 653 "doi = {10.1103/PhysRevC.82.024313}\n}\n",&set);CHKERRQ(ierr); 654 tao->niter=0; 655 if (tao->XL && tao->XU) { 656 /* Check x0 <= XU */ 657 PetscReal val; 658 659 CHKERRQ(VecCopy(tao->solution,mfqP->Xhist[0])); 660 CHKERRQ(VecAXPY(mfqP->Xhist[0],-1.0,tao->XU)); 661 CHKERRQ(VecMax(mfqP->Xhist[0],NULL,&val)); 662 PetscCheck(val <= 1e-10,PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_OUTOFRANGE,"X0 > upper bound"); 663 664 /* Check x0 >= xl */ 665 CHKERRQ(VecCopy(tao->XL,mfqP->Xhist[0])); 666 CHKERRQ(VecAXPY(mfqP->Xhist[0],-1.0,tao->solution)); 667 CHKERRQ(VecMax(mfqP->Xhist[0],NULL,&val)); 668 PetscCheck(val <= 1e-10,PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_OUTOFRANGE,"X0 < lower bound"); 669 670 /* Check x0 + delta < XU -- should be able to get around this eventually */ 671 672 CHKERRQ(VecSet(mfqP->Xhist[0],mfqP->delta)); 673 CHKERRQ(VecAXPY(mfqP->Xhist[0],1.0,tao->solution)); 674 CHKERRQ(VecAXPY(mfqP->Xhist[0],-1.0,tao->XU)); 675 CHKERRQ(VecMax(mfqP->Xhist[0],NULL,&val)); 676 PetscCheck(val <= 1e-10,PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_OUTOFRANGE,"X0 + delta > upper bound"); 677 } 678 679 blasm = mfqP->m; blasn=mfqP->n; blasnpmax = mfqP->npmax; 680 for (i=0;i<mfqP->n*mfqP->n*mfqP->m;++i) mfqP->H[i]=0; 681 682 CHKERRQ(VecCopy(tao->solution,mfqP->Xhist[0])); 683 684 /* This provides enough information to approximate the gradient of the objective */ 685 /* using a forward difference scheme. */ 686 687 CHKERRQ(PetscInfo(tao,"Initialize simplex; delta = %10.9e\n",(double)mfqP->delta)); 688 CHKERRQ(pounders_feval(tao,mfqP->Xhist[0],mfqP->Fhist[0],&mfqP->Fres[0])); 689 mfqP->minindex = 0; 690 minnorm = mfqP->Fres[0]; 691 692 CHKERRQ(VecGetOwnershipRange(mfqP->Xhist[0],&low,&high)); 693 for (i=1;i<mfqP->n+1;++i) { 694 CHKERRQ(VecCopy(mfqP->Xhist[0],mfqP->Xhist[i])); 695 696 if (i-1 >= low && i-1 < high) { 697 CHKERRQ(VecGetArray(mfqP->Xhist[i],&x)); 698 x[i-1-low] += mfqP->delta; 699 CHKERRQ(VecRestoreArray(mfqP->Xhist[i],&x)); 700 } 701 CHKMEMQ; 702 CHKERRQ(pounders_feval(tao,mfqP->Xhist[i],mfqP->Fhist[i],&mfqP->Fres[i])); 703 if (mfqP->Fres[i] < minnorm) { 704 mfqP->minindex = i; 705 minnorm = mfqP->Fres[i]; 706 } 707 } 708 CHKERRQ(VecCopy(mfqP->Xhist[mfqP->minindex],tao->solution)); 709 CHKERRQ(VecCopy(mfqP->Fhist[mfqP->minindex],tao->ls_res)); 710 CHKERRQ(PetscInfo(tao,"Finalize simplex; minnorm = %10.9e\n",(double)minnorm)); 711 712 /* Gather mpi vecs to one big local vec */ 713 714 /* Begin serial code */ 715 716 /* Disp[i] = Xi-xmin, i=1,..,mfqP->minindex-1,mfqP->minindex+1,..,n */ 717 /* Fdiff[i] = (Fi-Fmin)', i=1,..,mfqP->minindex-1,mfqP->minindex+1,..,n */ 718 /* (Column oriented for blas calls) */ 719 ii=0; 720 721 CHKERRQ(PetscInfo(tao,"Build matrix: %D\n",(PetscInt)mfqP->size)); 722 if (1 == mfqP->size) { 723 CHKERRQ(VecGetArrayRead(mfqP->Xhist[mfqP->minindex],&xmint)); 724 for (i=0;i<mfqP->n;i++) mfqP->xmin[i] = xmint[i]; 725 CHKERRQ(VecRestoreArrayRead(mfqP->Xhist[mfqP->minindex],&xmint)); 726 CHKERRQ(VecGetArrayRead(mfqP->Fhist[mfqP->minindex],&fmin)); 727 for (i=0;i<mfqP->n+1;i++) { 728 if (i == mfqP->minindex) continue; 729 730 CHKERRQ(VecGetArray(mfqP->Xhist[i],&x)); 731 for (j=0;j<mfqP->n;j++) { 732 mfqP->Disp[ii+mfqP->npmax*j] = (x[j] - mfqP->xmin[j])/mfqP->delta; 733 } 734 CHKERRQ(VecRestoreArray(mfqP->Xhist[i],&x)); 735 736 CHKERRQ(VecGetArray(mfqP->Fhist[i],&f)); 737 for (j=0;j<mfqP->m;j++) { 738 mfqP->Fdiff[ii+mfqP->n*j] = f[j] - fmin[j]; 739 } 740 CHKERRQ(VecRestoreArray(mfqP->Fhist[i],&f)); 741 742 mfqP->model_indices[ii++] = i; 743 } 744 for (j=0;j<mfqP->m;j++) { 745 mfqP->C[j] = fmin[j]; 746 } 747 CHKERRQ(VecRestoreArrayRead(mfqP->Fhist[mfqP->minindex],&fmin)); 748 } else { 749 CHKERRQ(VecSet(mfqP->localxmin,0)); 750 CHKERRQ(VecScatterBegin(mfqP->scatterx,mfqP->Xhist[mfqP->minindex],mfqP->localxmin,INSERT_VALUES,SCATTER_FORWARD)); 751 CHKERRQ(VecScatterEnd(mfqP->scatterx,mfqP->Xhist[mfqP->minindex],mfqP->localxmin,INSERT_VALUES,SCATTER_FORWARD)); 752 753 CHKERRQ(VecGetArrayRead(mfqP->localxmin,&xmint)); 754 for (i=0;i<mfqP->n;i++) mfqP->xmin[i] = xmint[i]; 755 CHKERRQ(VecRestoreArrayRead(mfqP->localxmin,&xmint)); 756 757 CHKERRQ(VecScatterBegin(mfqP->scatterf,mfqP->Fhist[mfqP->minindex],mfqP->localfmin,INSERT_VALUES,SCATTER_FORWARD)); 758 CHKERRQ(VecScatterEnd(mfqP->scatterf,mfqP->Fhist[mfqP->minindex],mfqP->localfmin,INSERT_VALUES,SCATTER_FORWARD)); 759 CHKERRQ(VecGetArrayRead(mfqP->localfmin,&fmin)); 760 for (i=0;i<mfqP->n+1;i++) { 761 if (i == mfqP->minindex) continue; 762 763 CHKERRQ(VecScatterBegin(mfqP->scatterx,mfqP->Xhist[ii],mfqP->localx,INSERT_VALUES, SCATTER_FORWARD)); 764 CHKERRQ(VecScatterEnd(mfqP->scatterx,mfqP->Xhist[ii],mfqP->localx,INSERT_VALUES, SCATTER_FORWARD)); 765 CHKERRQ(VecGetArray(mfqP->localx,&x)); 766 for (j=0;j<mfqP->n;j++) { 767 mfqP->Disp[ii+mfqP->npmax*j] = (x[j] - mfqP->xmin[j])/mfqP->delta; 768 } 769 CHKERRQ(VecRestoreArray(mfqP->localx,&x)); 770 771 CHKERRQ(VecScatterBegin(mfqP->scatterf,mfqP->Fhist[ii],mfqP->localf,INSERT_VALUES, SCATTER_FORWARD)); 772 CHKERRQ(VecScatterEnd(mfqP->scatterf,mfqP->Fhist[ii],mfqP->localf,INSERT_VALUES, SCATTER_FORWARD)); 773 CHKERRQ(VecGetArray(mfqP->localf,&f)); 774 for (j=0;j<mfqP->m;j++) { 775 mfqP->Fdiff[ii+mfqP->n*j] = f[j] - fmin[j]; 776 } 777 CHKERRQ(VecRestoreArray(mfqP->localf,&f)); 778 779 mfqP->model_indices[ii++] = i; 780 } 781 for (j=0;j<mfqP->m;j++) { 782 mfqP->C[j] = fmin[j]; 783 } 784 CHKERRQ(VecRestoreArrayRead(mfqP->localfmin,&fmin)); 785 } 786 787 /* Determine the initial quadratic models */ 788 /* G = D(ModelIn,:) \ (F(ModelIn,1:m)-repmat(F(xkin,1:m),n,1)); */ 789 /* D (nxn) Fdiff (nxm) => G (nxm) */ 790 blasncopy = blasn; 791 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&blasn,&blasm,mfqP->Disp,&blasnpmax,mfqP->iwork,mfqP->Fdiff,&blasncopy,&info)); 792 CHKERRQ(PetscInfo(tao,"Linear solve return: %D\n",(PetscInt)info)); 793 794 CHKERRQ(pounders_update_res(tao)); 795 796 valid = PETSC_TRUE; 797 798 CHKERRQ(VecSetValues(tao->gradient,mfqP->n,mfqP->indices,mfqP->Gres,INSERT_VALUES)); 799 CHKERRQ(VecAssemblyBegin(tao->gradient)); 800 CHKERRQ(VecAssemblyEnd(tao->gradient)); 801 CHKERRQ(VecNorm(tao->gradient,NORM_2,&gnorm)); 802 gnorm *= mfqP->delta; 803 CHKERRQ(VecCopy(mfqP->Xhist[mfqP->minindex],tao->solution)); 804 805 tao->reason = TAO_CONTINUE_ITERATING; 806 CHKERRQ(TaoLogConvergenceHistory(tao,minnorm,gnorm,0.0,tao->ksp_its)); 807 CHKERRQ(TaoMonitor(tao,tao->niter,minnorm,gnorm,0.0,step)); 808 CHKERRQ((*tao->ops->convergencetest)(tao,tao->cnvP)); 809 810 mfqP->nHist = mfqP->n+1; 811 mfqP->nmodelpoints = mfqP->n+1; 812 CHKERRQ(PetscInfo(tao,"Initial gradient: %20.19e\n",(double)gnorm)); 813 814 while (tao->reason == TAO_CONTINUE_ITERATING) { 815 PetscReal gnm = 1e-4; 816 /* Call general purpose update function */ 817 if (tao->ops->update) { 818 CHKERRQ((*tao->ops->update)(tao, tao->niter, tao->user_update)); 819 } 820 tao->niter++; 821 /* Solve the subproblem min{Q(s): ||s|| <= 1.0} */ 822 CHKERRQ(gqtwrap(tao,&gnm,&mdec)); 823 /* Evaluate the function at the new point */ 824 825 for (i=0;i<mfqP->n;i++) { 826 mfqP->work[i] = mfqP->Xsubproblem[i]*mfqP->delta + mfqP->xmin[i]; 827 } 828 CHKERRQ(VecDuplicate(tao->solution,&mfqP->Xhist[mfqP->nHist])); 829 CHKERRQ(VecDuplicate(tao->ls_res,&mfqP->Fhist[mfqP->nHist])); 830 CHKERRQ(VecSetValues(mfqP->Xhist[mfqP->nHist],mfqP->n,mfqP->indices,mfqP->work,INSERT_VALUES)); 831 CHKERRQ(VecAssemblyBegin(mfqP->Xhist[mfqP->nHist])); 832 CHKERRQ(VecAssemblyEnd(mfqP->Xhist[mfqP->nHist])); 833 834 CHKERRQ(pounders_feval(tao,mfqP->Xhist[mfqP->nHist],mfqP->Fhist[mfqP->nHist],&mfqP->Fres[mfqP->nHist])); 835 836 rho = (mfqP->Fres[mfqP->minindex] - mfqP->Fres[mfqP->nHist]) / mdec; 837 mfqP->nHist++; 838 839 /* Update the center */ 840 if ((rho >= mfqP->eta1) || (rho > mfqP->eta0 && valid==PETSC_TRUE)) { 841 /* Update model to reflect new base point */ 842 for (i=0;i<mfqP->n;i++) { 843 mfqP->work[i] = (mfqP->work[i] - mfqP->xmin[i])/mfqP->delta; 844 } 845 for (j=0;j<mfqP->m;j++) { 846 /* C(j) = C(j) + work*G(:,j) + .5*work*H(:,:,j)*work'; 847 G(:,j) = G(:,j) + H(:,:,j)*work' */ 848 for (k=0;k<mfqP->n;k++) { 849 mfqP->work2[k]=0.0; 850 for (l=0;l<mfqP->n;l++) { 851 mfqP->work2[k]+=mfqP->H[j + mfqP->m*(k + l*mfqP->n)]*mfqP->work[l]; 852 } 853 } 854 for (i=0;i<mfqP->n;i++) { 855 mfqP->C[j]+=mfqP->work[i]*(mfqP->Fdiff[i + mfqP->n* j] + 0.5*mfqP->work2[i]); 856 mfqP->Fdiff[i+mfqP->n*j] +=mfqP-> work2[i]; 857 } 858 } 859 /* Cres += work*Gres + .5*work*Hres*work'; 860 Gres += Hres*work'; */ 861 862 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&blasn,&blasn,&one,mfqP->Hres,&blasn,mfqP->work,&ione,&zero,mfqP->work2,&ione)); 863 for (i=0;i<mfqP->n;i++) { 864 mfqP->Gres[i] += mfqP->work2[i]; 865 } 866 mfqP->minindex = mfqP->nHist-1; 867 minnorm = mfqP->Fres[mfqP->minindex]; 868 CHKERRQ(VecCopy(mfqP->Fhist[mfqP->minindex],tao->ls_res)); 869 /* Change current center */ 870 CHKERRQ(VecGetArrayRead(mfqP->Xhist[mfqP->minindex],&xmint)); 871 for (i=0;i<mfqP->n;i++) { 872 mfqP->xmin[i] = xmint[i]; 873 } 874 CHKERRQ(VecRestoreArrayRead(mfqP->Xhist[mfqP->minindex],&xmint)); 875 } 876 877 /* Evaluate at a model-improving point if necessary */ 878 if (valid == PETSC_FALSE) { 879 mfqP->q_is_I = 1; 880 mfqP->nmodelpoints = 0; 881 CHKERRQ(affpoints(mfqP,mfqP->xmin,mfqP->c1)); 882 if (mfqP->nmodelpoints < mfqP->n) { 883 CHKERRQ(PetscInfo(tao,"Model not valid -- model-improving\n")); 884 CHKERRQ(modelimprove(tao,mfqP,1)); 885 } 886 } 887 888 /* Update the trust region radius */ 889 deltaold = mfqP->delta; 890 normxsp = 0; 891 for (i=0;i<mfqP->n;i++) { 892 normxsp += mfqP->Xsubproblem[i]*mfqP->Xsubproblem[i]; 893 } 894 normxsp = PetscSqrtReal(normxsp); 895 if (rho >= mfqP->eta1 && normxsp > 0.5*mfqP->delta) { 896 mfqP->delta = PetscMin(mfqP->delta*mfqP->gamma1,mfqP->deltamax); 897 } else if (valid == PETSC_TRUE) { 898 mfqP->delta = PetscMax(mfqP->delta*mfqP->gamma0,mfqP->deltamin); 899 } 900 901 /* Compute the next interpolation set */ 902 mfqP->q_is_I = 1; 903 mfqP->nmodelpoints=0; 904 CHKERRQ(PetscInfo(tao,"Affine Points: xmin = %20.19e, c1 = %20.19e\n",(double)*mfqP->xmin,(double)mfqP->c1)); 905 CHKERRQ(affpoints(mfqP,mfqP->xmin,mfqP->c1)); 906 if (mfqP->nmodelpoints == mfqP->n) { 907 valid = PETSC_TRUE; 908 } else { 909 valid = PETSC_FALSE; 910 CHKERRQ(PetscInfo(tao,"Affine Points: xmin = %20.19e, c2 = %20.19e\n",(double)*mfqP->xmin,(double)mfqP->c2)); 911 CHKERRQ(affpoints(mfqP,mfqP->xmin,mfqP->c2)); 912 if (mfqP->n > mfqP->nmodelpoints) { 913 CHKERRQ(PetscInfo(tao,"Model not valid -- adding geometry points\n")); 914 CHKERRQ(modelimprove(tao,mfqP,mfqP->n - mfqP->nmodelpoints)); 915 } 916 } 917 for (i=mfqP->nmodelpoints;i>0;i--) { 918 mfqP->model_indices[i] = mfqP->model_indices[i-1]; 919 } 920 mfqP->nmodelpoints++; 921 mfqP->model_indices[0] = mfqP->minindex; 922 CHKERRQ(morepoints(mfqP)); 923 for (i=0;i<mfqP->nmodelpoints;i++) { 924 CHKERRQ(VecGetArray(mfqP->Xhist[mfqP->model_indices[i]],&x)); 925 for (j=0;j<mfqP->n;j++) { 926 mfqP->Disp[i + mfqP->npmax*j] = (x[j] - mfqP->xmin[j]) / deltaold; 927 } 928 CHKERRQ(VecRestoreArray(mfqP->Xhist[mfqP->model_indices[i]],&x)); 929 CHKERRQ(VecGetArray(mfqP->Fhist[mfqP->model_indices[i]],&f)); 930 for (j=0;j<mfqP->m;j++) { 931 for (k=0;k<mfqP->n;k++) { 932 mfqP->work[k]=0.0; 933 for (l=0;l<mfqP->n;l++) { 934 mfqP->work[k] += mfqP->H[j + mfqP->m*(k + mfqP->n*l)] * mfqP->Disp[i + mfqP->npmax*l]; 935 } 936 } 937 PetscStackCallBLAS("BLASdot",mfqP->RES[j*mfqP->npmax + i] = -mfqP->C[j] - BLASdot_(&blasn,&mfqP->Fdiff[j*mfqP->n],&ione,&mfqP->Disp[i],&blasnpmax) - 0.5*BLASdot_(&blasn,mfqP->work,&ione,&mfqP->Disp[i],&blasnpmax) + f[j]); 938 } 939 CHKERRQ(VecRestoreArray(mfqP->Fhist[mfqP->model_indices[i]],&f)); 940 } 941 942 /* Update the quadratic model */ 943 CHKERRQ(PetscInfo(tao,"Get Quad, size: %D, points: %D\n",mfqP->n,mfqP->nmodelpoints)); 944 CHKERRQ(getquadpounders(mfqP)); 945 CHKERRQ(VecGetArrayRead(mfqP->Fhist[mfqP->minindex],&fmin)); 946 PetscStackCallBLAS("BLAScopy",BLAScopy_(&blasm,fmin,&ione,mfqP->C,&ione)); 947 /* G = G*(delta/deltaold) + Gdel */ 948 ratio = mfqP->delta/deltaold; 949 iblas = blasm*blasn; 950 PetscStackCallBLAS("BLASscal",BLASscal_(&iblas,&ratio,mfqP->Fdiff,&ione)); 951 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&iblas,&one,mfqP->Gdel,&ione,mfqP->Fdiff,&ione)); 952 /* H = H*(delta/deltaold)^2 + Hdel */ 953 iblas = blasm*blasn*blasn; 954 ratio *= ratio; 955 PetscStackCallBLAS("BLASscal",BLASscal_(&iblas,&ratio,mfqP->H,&ione)); 956 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&iblas,&one,mfqP->Hdel,&ione,mfqP->H,&ione)); 957 958 /* Get residuals */ 959 CHKERRQ(pounders_update_res(tao)); 960 961 /* Export solution and gradient residual to TAO */ 962 CHKERRQ(VecCopy(mfqP->Xhist[mfqP->minindex],tao->solution)); 963 CHKERRQ(VecSetValues(tao->gradient,mfqP->n,mfqP->indices,mfqP->Gres,INSERT_VALUES)); 964 CHKERRQ(VecAssemblyBegin(tao->gradient)); 965 CHKERRQ(VecAssemblyEnd(tao->gradient)); 966 CHKERRQ(VecNorm(tao->gradient,NORM_2,&gnorm)); 967 gnorm *= mfqP->delta; 968 /* final criticality test */ 969 CHKERRQ(TaoLogConvergenceHistory(tao,minnorm,gnorm,0.0,tao->ksp_its)); 970 CHKERRQ(TaoMonitor(tao,tao->niter,minnorm,gnorm,0.0,step)); 971 CHKERRQ((*tao->ops->convergencetest)(tao,tao->cnvP)); 972 /* test for repeated model */ 973 if (mfqP->nmodelpoints==mfqP->last_nmodelpoints) { 974 same = PETSC_TRUE; 975 } else { 976 same = PETSC_FALSE; 977 } 978 for (i=0;i<mfqP->nmodelpoints;i++) { 979 if (same) { 980 if (mfqP->model_indices[i] == mfqP->last_model_indices[i]) { 981 same = PETSC_TRUE; 982 } else { 983 same = PETSC_FALSE; 984 } 985 } 986 mfqP->last_model_indices[i] = mfqP->model_indices[i]; 987 } 988 mfqP->last_nmodelpoints = mfqP->nmodelpoints; 989 if (same && mfqP->delta == deltaold) { 990 CHKERRQ(PetscInfo(tao,"Identical model used in successive iterations\n")); 991 tao->reason = TAO_CONVERGED_STEPTOL; 992 } 993 } 994 PetscFunctionReturn(0); 995 } 996 997 static PetscErrorCode TaoSetUp_POUNDERS(Tao tao) 998 { 999 TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data; 1000 PetscInt i,j; 1001 IS isfloc,isfglob,isxloc,isxglob; 1002 1003 PetscFunctionBegin; 1004 if (!tao->gradient) CHKERRQ(VecDuplicate(tao->solution,&tao->gradient)); 1005 if (!tao->stepdirection) CHKERRQ(VecDuplicate(tao->solution,&tao->stepdirection)); 1006 CHKERRQ(VecGetSize(tao->solution,&mfqP->n)); 1007 CHKERRQ(VecGetSize(tao->ls_res,&mfqP->m)); 1008 mfqP->c1 = PetscSqrtReal((PetscReal)mfqP->n); 1009 if (mfqP->npmax == PETSC_DEFAULT) { 1010 mfqP->npmax = 2*mfqP->n + 1; 1011 } 1012 mfqP->npmax = PetscMin((mfqP->n+1)*(mfqP->n+2)/2,mfqP->npmax); 1013 mfqP->npmax = PetscMax(mfqP->npmax, mfqP->n+2); 1014 1015 CHKERRQ(PetscMalloc1(tao->max_funcs+100,&mfqP->Xhist)); 1016 CHKERRQ(PetscMalloc1(tao->max_funcs+100,&mfqP->Fhist)); 1017 for (i=0;i<mfqP->n+1;i++) { 1018 CHKERRQ(VecDuplicate(tao->solution,&mfqP->Xhist[i])); 1019 CHKERRQ(VecDuplicate(tao->ls_res,&mfqP->Fhist[i])); 1020 } 1021 CHKERRQ(VecDuplicate(tao->solution,&mfqP->workxvec)); 1022 CHKERRQ(VecDuplicate(tao->ls_res,&mfqP->workfvec)); 1023 mfqP->nHist = 0; 1024 1025 CHKERRQ(PetscMalloc1(tao->max_funcs+100,&mfqP->Fres)); 1026 CHKERRQ(PetscMalloc1(mfqP->npmax*mfqP->m,&mfqP->RES)); 1027 CHKERRQ(PetscMalloc1(mfqP->n,&mfqP->work)); 1028 CHKERRQ(PetscMalloc1(mfqP->n,&mfqP->work2)); 1029 CHKERRQ(PetscMalloc1(mfqP->n,&mfqP->work3)); 1030 CHKERRQ(PetscMalloc1(PetscMax(mfqP->m,mfqP->n+1),&mfqP->mwork)); 1031 CHKERRQ(PetscMalloc1(mfqP->npmax - mfqP->n - 1,&mfqP->omega)); 1032 CHKERRQ(PetscMalloc1(mfqP->n * (mfqP->n+1) / 2,&mfqP->beta)); 1033 CHKERRQ(PetscMalloc1(mfqP->n + 1 ,&mfqP->alpha)); 1034 1035 CHKERRQ(PetscMalloc1(mfqP->n*mfqP->n*mfqP->m,&mfqP->H)); 1036 CHKERRQ(PetscMalloc1(mfqP->npmax*mfqP->npmax,&mfqP->Q)); 1037 CHKERRQ(PetscMalloc1(mfqP->npmax*mfqP->npmax,&mfqP->Q_tmp)); 1038 CHKERRQ(PetscMalloc1(mfqP->n*(mfqP->n+1)/2*(mfqP->npmax),&mfqP->L)); 1039 CHKERRQ(PetscMalloc1(mfqP->n*(mfqP->n+1)/2*(mfqP->npmax),&mfqP->L_tmp)); 1040 CHKERRQ(PetscMalloc1(mfqP->n*(mfqP->n+1)/2*(mfqP->npmax),&mfqP->L_save)); 1041 CHKERRQ(PetscMalloc1(mfqP->n*(mfqP->n+1)/2*(mfqP->npmax),&mfqP->N)); 1042 CHKERRQ(PetscMalloc1(mfqP->npmax * (mfqP->n+1) ,&mfqP->M)); 1043 CHKERRQ(PetscMalloc1(mfqP->npmax * (mfqP->npmax - mfqP->n - 1) , &mfqP->Z)); 1044 CHKERRQ(PetscMalloc1(mfqP->npmax,&mfqP->tau)); 1045 CHKERRQ(PetscMalloc1(mfqP->npmax,&mfqP->tau_tmp)); 1046 mfqP->nmax = PetscMax(5*mfqP->npmax,mfqP->n*(mfqP->n+1)/2); 1047 CHKERRQ(PetscMalloc1(mfqP->nmax,&mfqP->npmaxwork)); 1048 CHKERRQ(PetscMalloc1(mfqP->nmax,&mfqP->npmaxiwork)); 1049 CHKERRQ(PetscMalloc1(mfqP->n,&mfqP->xmin)); 1050 CHKERRQ(PetscMalloc1(mfqP->m,&mfqP->C)); 1051 CHKERRQ(PetscMalloc1(mfqP->m*mfqP->n,&mfqP->Fdiff)); 1052 CHKERRQ(PetscMalloc1(mfqP->npmax*mfqP->n,&mfqP->Disp)); 1053 CHKERRQ(PetscMalloc1(mfqP->n,&mfqP->Gres)); 1054 CHKERRQ(PetscMalloc1(mfqP->n*mfqP->n,&mfqP->Hres)); 1055 CHKERRQ(PetscMalloc1(mfqP->n*mfqP->n,&mfqP->Gpoints)); 1056 CHKERRQ(PetscMalloc1(mfqP->npmax,&mfqP->model_indices)); 1057 CHKERRQ(PetscMalloc1(mfqP->npmax,&mfqP->last_model_indices)); 1058 CHKERRQ(PetscMalloc1(mfqP->n,&mfqP->Xsubproblem)); 1059 CHKERRQ(PetscMalloc1(mfqP->m*mfqP->n,&mfqP->Gdel)); 1060 CHKERRQ(PetscMalloc1(mfqP->n*mfqP->n*mfqP->m, &mfqP->Hdel)); 1061 CHKERRQ(PetscMalloc1(PetscMax(mfqP->m,mfqP->n),&mfqP->indices)); 1062 CHKERRQ(PetscMalloc1(mfqP->n,&mfqP->iwork)); 1063 CHKERRQ(PetscMalloc1(mfqP->m*mfqP->m,&mfqP->w)); 1064 for (i=0;i<mfqP->m;i++) { 1065 for (j=0;j<mfqP->m;j++) { 1066 if (i==j) { 1067 mfqP->w[i+mfqP->m*j]=1.0; 1068 } else { 1069 mfqP->w[i+mfqP->m*j]=0.0; 1070 } 1071 } 1072 } 1073 for (i=0;i<PetscMax(mfqP->m,mfqP->n);i++) { 1074 mfqP->indices[i] = i; 1075 } 1076 CHKERRMPI(MPI_Comm_size(((PetscObject)tao)->comm,&mfqP->size)); 1077 if (mfqP->size > 1) { 1078 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,mfqP->n,&mfqP->localx)); 1079 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,mfqP->n,&mfqP->localxmin)); 1080 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,mfqP->m,&mfqP->localf)); 1081 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,mfqP->m,&mfqP->localfmin)); 1082 CHKERRQ(ISCreateStride(MPI_COMM_SELF,mfqP->n,0,1,&isxloc)); 1083 CHKERRQ(ISCreateStride(MPI_COMM_SELF,mfqP->n,0,1,&isxglob)); 1084 CHKERRQ(ISCreateStride(MPI_COMM_SELF,mfqP->m,0,1,&isfloc)); 1085 CHKERRQ(ISCreateStride(MPI_COMM_SELF,mfqP->m,0,1,&isfglob)); 1086 1087 CHKERRQ(VecScatterCreate(tao->solution,isxglob,mfqP->localx,isxloc,&mfqP->scatterx)); 1088 CHKERRQ(VecScatterCreate(tao->ls_res,isfglob,mfqP->localf,isfloc,&mfqP->scatterf)); 1089 1090 CHKERRQ(ISDestroy(&isxloc)); 1091 CHKERRQ(ISDestroy(&isxglob)); 1092 CHKERRQ(ISDestroy(&isfloc)); 1093 CHKERRQ(ISDestroy(&isfglob)); 1094 } 1095 1096 if (!mfqP->usegqt) { 1097 KSP ksp; 1098 PC pc; 1099 CHKERRQ(VecCreateSeqWithArray(PETSC_COMM_SELF,mfqP->n,mfqP->n,mfqP->Xsubproblem,&mfqP->subx)); 1100 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,mfqP->n,&mfqP->subxl)); 1101 CHKERRQ(VecDuplicate(mfqP->subxl,&mfqP->subb)); 1102 CHKERRQ(VecDuplicate(mfqP->subxl,&mfqP->subxu)); 1103 CHKERRQ(VecDuplicate(mfqP->subxl,&mfqP->subpdel)); 1104 CHKERRQ(VecDuplicate(mfqP->subxl,&mfqP->subndel)); 1105 CHKERRQ(TaoCreate(PETSC_COMM_SELF,&mfqP->subtao)); 1106 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)mfqP->subtao, (PetscObject)tao, 1)); 1107 CHKERRQ(TaoSetType(mfqP->subtao,TAOBNTR)); 1108 CHKERRQ(TaoSetOptionsPrefix(mfqP->subtao,"pounders_subsolver_")); 1109 CHKERRQ(TaoSetSolution(mfqP->subtao,mfqP->subx)); 1110 CHKERRQ(TaoSetObjectiveAndGradient(mfqP->subtao,NULL,pounders_fg,(void*)mfqP)); 1111 CHKERRQ(TaoSetMaximumIterations(mfqP->subtao,mfqP->gqt_maxits)); 1112 CHKERRQ(TaoSetFromOptions(mfqP->subtao)); 1113 CHKERRQ(TaoGetKSP(mfqP->subtao,&ksp)); 1114 if (ksp) { 1115 CHKERRQ(KSPGetPC(ksp,&pc)); 1116 CHKERRQ(PCSetType(pc,PCNONE)); 1117 } 1118 CHKERRQ(TaoSetVariableBounds(mfqP->subtao,mfqP->subxl,mfqP->subxu)); 1119 CHKERRQ(MatCreateSeqDense(PETSC_COMM_SELF,mfqP->n,mfqP->n,mfqP->Hres,&mfqP->subH)); 1120 CHKERRQ(TaoSetHessian(mfqP->subtao,mfqP->subH,mfqP->subH,pounders_h,(void*)mfqP)); 1121 } 1122 PetscFunctionReturn(0); 1123 } 1124 1125 static PetscErrorCode TaoDestroy_POUNDERS(Tao tao) 1126 { 1127 TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data; 1128 PetscInt i; 1129 1130 PetscFunctionBegin; 1131 if (!mfqP->usegqt) { 1132 CHKERRQ(TaoDestroy(&mfqP->subtao)); 1133 CHKERRQ(VecDestroy(&mfqP->subx)); 1134 CHKERRQ(VecDestroy(&mfqP->subxl)); 1135 CHKERRQ(VecDestroy(&mfqP->subxu)); 1136 CHKERRQ(VecDestroy(&mfqP->subb)); 1137 CHKERRQ(VecDestroy(&mfqP->subpdel)); 1138 CHKERRQ(VecDestroy(&mfqP->subndel)); 1139 CHKERRQ(MatDestroy(&mfqP->subH)); 1140 } 1141 CHKERRQ(PetscFree(mfqP->Fres)); 1142 CHKERRQ(PetscFree(mfqP->RES)); 1143 CHKERRQ(PetscFree(mfqP->work)); 1144 CHKERRQ(PetscFree(mfqP->work2)); 1145 CHKERRQ(PetscFree(mfqP->work3)); 1146 CHKERRQ(PetscFree(mfqP->mwork)); 1147 CHKERRQ(PetscFree(mfqP->omega)); 1148 CHKERRQ(PetscFree(mfqP->beta)); 1149 CHKERRQ(PetscFree(mfqP->alpha)); 1150 CHKERRQ(PetscFree(mfqP->H)); 1151 CHKERRQ(PetscFree(mfqP->Q)); 1152 CHKERRQ(PetscFree(mfqP->Q_tmp)); 1153 CHKERRQ(PetscFree(mfqP->L)); 1154 CHKERRQ(PetscFree(mfqP->L_tmp)); 1155 CHKERRQ(PetscFree(mfqP->L_save)); 1156 CHKERRQ(PetscFree(mfqP->N)); 1157 CHKERRQ(PetscFree(mfqP->M)); 1158 CHKERRQ(PetscFree(mfqP->Z)); 1159 CHKERRQ(PetscFree(mfqP->tau)); 1160 CHKERRQ(PetscFree(mfqP->tau_tmp)); 1161 CHKERRQ(PetscFree(mfqP->npmaxwork)); 1162 CHKERRQ(PetscFree(mfqP->npmaxiwork)); 1163 CHKERRQ(PetscFree(mfqP->xmin)); 1164 CHKERRQ(PetscFree(mfqP->C)); 1165 CHKERRQ(PetscFree(mfqP->Fdiff)); 1166 CHKERRQ(PetscFree(mfqP->Disp)); 1167 CHKERRQ(PetscFree(mfqP->Gres)); 1168 CHKERRQ(PetscFree(mfqP->Hres)); 1169 CHKERRQ(PetscFree(mfqP->Gpoints)); 1170 CHKERRQ(PetscFree(mfqP->model_indices)); 1171 CHKERRQ(PetscFree(mfqP->last_model_indices)); 1172 CHKERRQ(PetscFree(mfqP->Xsubproblem)); 1173 CHKERRQ(PetscFree(mfqP->Gdel)); 1174 CHKERRQ(PetscFree(mfqP->Hdel)); 1175 CHKERRQ(PetscFree(mfqP->indices)); 1176 CHKERRQ(PetscFree(mfqP->iwork)); 1177 CHKERRQ(PetscFree(mfqP->w)); 1178 for (i=0;i<mfqP->nHist;i++) { 1179 CHKERRQ(VecDestroy(&mfqP->Xhist[i])); 1180 CHKERRQ(VecDestroy(&mfqP->Fhist[i])); 1181 } 1182 CHKERRQ(VecDestroy(&mfqP->workxvec)); 1183 CHKERRQ(VecDestroy(&mfqP->workfvec)); 1184 CHKERRQ(PetscFree(mfqP->Xhist)); 1185 CHKERRQ(PetscFree(mfqP->Fhist)); 1186 1187 if (mfqP->size > 1) { 1188 CHKERRQ(VecDestroy(&mfqP->localx)); 1189 CHKERRQ(VecDestroy(&mfqP->localxmin)); 1190 CHKERRQ(VecDestroy(&mfqP->localf)); 1191 CHKERRQ(VecDestroy(&mfqP->localfmin)); 1192 } 1193 CHKERRQ(PetscFree(tao->data)); 1194 PetscFunctionReturn(0); 1195 } 1196 1197 static PetscErrorCode TaoSetFromOptions_POUNDERS(PetscOptionItems *PetscOptionsObject,Tao tao) 1198 { 1199 TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data; 1200 1201 PetscFunctionBegin; 1202 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"POUNDERS method for least-squares optimization")); 1203 CHKERRQ(PetscOptionsReal("-tao_pounders_delta","initial delta","",mfqP->delta,&mfqP->delta0,NULL)); 1204 mfqP->delta = mfqP->delta0; 1205 CHKERRQ(PetscOptionsInt("-tao_pounders_npmax","max number of points in model","",mfqP->npmax,&mfqP->npmax,NULL)); 1206 CHKERRQ(PetscOptionsBool("-tao_pounders_gqt","use gqt algorithm for subproblem","",mfqP->usegqt,&mfqP->usegqt,NULL)); 1207 CHKERRQ(PetscOptionsTail()); 1208 PetscFunctionReturn(0); 1209 } 1210 1211 static PetscErrorCode TaoView_POUNDERS(Tao tao, PetscViewer viewer) 1212 { 1213 TAO_POUNDERS *mfqP = (TAO_POUNDERS *)tao->data; 1214 PetscBool isascii; 1215 1216 PetscFunctionBegin; 1217 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 1218 if (isascii) { 1219 CHKERRQ(PetscViewerASCIIPrintf(viewer, "initial delta: %g\n",(double)mfqP->delta0)); 1220 CHKERRQ(PetscViewerASCIIPrintf(viewer, "final delta: %g\n",(double)mfqP->delta)); 1221 CHKERRQ(PetscViewerASCIIPrintf(viewer, "model points: %D\n",mfqP->nmodelpoints)); 1222 if (mfqP->usegqt) { 1223 CHKERRQ(PetscViewerASCIIPrintf(viewer, "subproblem solver: gqt\n")); 1224 } else { 1225 CHKERRQ(TaoView(mfqP->subtao, viewer)); 1226 } 1227 } 1228 PetscFunctionReturn(0); 1229 } 1230 /*MC 1231 TAOPOUNDERS - POUNDERS derivate-free model-based algorithm for nonlinear least squares 1232 1233 Options Database Keys: 1234 + -tao_pounders_delta - initial step length 1235 . -tao_pounders_npmax - maximum number of points in model 1236 - -tao_pounders_gqt - use gqt algorithm for subproblem instead of TRON 1237 1238 Level: beginner 1239 1240 M*/ 1241 1242 PETSC_EXTERN PetscErrorCode TaoCreate_POUNDERS(Tao tao) 1243 { 1244 TAO_POUNDERS *mfqP = (TAO_POUNDERS*)tao->data; 1245 1246 PetscFunctionBegin; 1247 tao->ops->setup = TaoSetUp_POUNDERS; 1248 tao->ops->solve = TaoSolve_POUNDERS; 1249 tao->ops->view = TaoView_POUNDERS; 1250 tao->ops->setfromoptions = TaoSetFromOptions_POUNDERS; 1251 tao->ops->destroy = TaoDestroy_POUNDERS; 1252 1253 CHKERRQ(PetscNewLog(tao,&mfqP)); 1254 tao->data = (void*)mfqP; 1255 /* Override default settings (unless already changed) */ 1256 if (!tao->max_it_changed) tao->max_it = 2000; 1257 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 1258 mfqP->npmax = PETSC_DEFAULT; 1259 mfqP->delta0 = 0.1; 1260 mfqP->delta = 0.1; 1261 mfqP->deltamax=1e3; 1262 mfqP->deltamin=1e-6; 1263 mfqP->c2 = 10.0; 1264 mfqP->theta1=1e-5; 1265 mfqP->theta2=1e-4; 1266 mfqP->gamma0=0.5; 1267 mfqP->gamma1=2.0; 1268 mfqP->eta0 = 0.0; 1269 mfqP->eta1 = 0.1; 1270 mfqP->usegqt = PETSC_FALSE; 1271 mfqP->gqt_rtol = 0.001; 1272 mfqP->gqt_maxits = 50; 1273 mfqP->workxvec = NULL; 1274 PetscFunctionReturn(0); 1275 } 1276