1 #include <petsctao.h> 2 /* 3 Description: ADMM tomography reconstruction example . 4 0.5*||Ax-b||^2 + lambda*g(x) 5 Reference: BRGN Tomography Example 6 */ 7 8 static char help[] = "Finds the ADMM solution to the under constraint linear model Ax = b, with regularizer. \n\ 9 A is a M*N real matrix (M<N), x is sparse. A good regularizer is an L1 regularizer. \n\ 10 We first split the operator into 0.5*||Ax-b||^2, f(x), and lambda*||x||_1, g(z), where lambda is user specified weight. \n\ 11 g(z) could be either ||z||_1, or ||z||_2^2. Default closed form solution for NORM1 would be soft-threshold, which is \n\ 12 natively supported in admm.c with -tao_admm_regularizer_type soft-threshold. Or user can use regular TAO solver for \n\ 13 either NORM1 or NORM2 or TAOSHELL, with -reg {1,2,3} \n\ 14 Then, we augment both f and g, and solve it via ADMM. \n\ 15 D is the M*N transform matrix so that D*x is sparse. \n"; 16 17 typedef struct { 18 PetscInt M,N,K,reg; 19 PetscReal lambda,eps,mumin; 20 Mat A,ATA,H,Hx,D,Hz,DTD,HF; 21 Vec c,xlb,xub,x,b,workM,workN,workN2,workN3,xGT; /* observation b, ground truth xGT, the lower bound and upper bound of x*/ 22 } AppCtx; 23 24 /*------------------------------------------------------------*/ 25 26 PetscErrorCode NullJacobian(Tao tao,Vec X,Mat J,Mat Jpre,void *ptr) 27 { 28 PetscFunctionBegin; 29 PetscFunctionReturn(0); 30 } 31 32 /*------------------------------------------------------------*/ 33 34 static PetscErrorCode TaoShellSolve_SoftThreshold(Tao tao) 35 { 36 PetscErrorCode ierr; 37 PetscReal lambda, mu; 38 AppCtx *user; 39 Vec out,work,y,x; 40 Tao admm_tao,misfit; 41 42 PetscFunctionBegin; 43 user = NULL; 44 mu = 0; 45 ierr = TaoGetADMMParentTao(tao,&admm_tao);CHKERRQ(ierr); 46 ierr = TaoADMMGetMisfitSubsolver(admm_tao, &misfit);CHKERRQ(ierr); 47 ierr = TaoADMMGetSpectralPenalty(admm_tao,&mu);CHKERRQ(ierr); 48 ierr = TaoShellGetContext(tao,&user);CHKERRQ(ierr); 49 50 lambda = user->lambda; 51 work = user->workN; 52 ierr = TaoGetSolutionVector(tao, &out);CHKERRQ(ierr); 53 ierr = TaoGetSolutionVector(misfit, &x);CHKERRQ(ierr); 54 ierr = TaoADMMGetDualVector(admm_tao, &y);CHKERRQ(ierr); 55 56 /* Dx + y/mu */ 57 ierr = MatMult(user->D,x,work);CHKERRQ(ierr); 58 ierr = VecAXPY(work,1/mu,y);CHKERRQ(ierr); 59 60 /* soft thresholding */ 61 ierr = TaoSoftThreshold(work, -lambda/mu, lambda/mu, out);CHKERRQ(ierr); 62 PetscFunctionReturn(0); 63 } 64 65 /*------------------------------------------------------------*/ 66 67 PetscErrorCode MisfitObjectiveAndGradient(Tao tao,Vec X,PetscReal *f,Vec g,void *ptr) 68 { 69 AppCtx *user = (AppCtx*)ptr; 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 /* Objective 0.5*||Ax-b||_2^2 */ 74 ierr = MatMult(user->A,X,user->workM);CHKERRQ(ierr); 75 ierr = VecAXPY(user->workM,-1,user->b);CHKERRQ(ierr); 76 ierr = VecDot(user->workM,user->workM,f);CHKERRQ(ierr); 77 *f *= 0.5; 78 /* Gradient. ATAx-ATb */ 79 ierr = MatMult(user->ATA,X,user->workN);CHKERRQ(ierr); 80 ierr = MatMultTranspose(user->A,user->b,user->workN2);CHKERRQ(ierr); 81 ierr = VecWAXPY(g,-1.,user->workN2,user->workN);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 /*------------------------------------------------------------*/ 86 87 PetscErrorCode RegularizerObjectiveAndGradient1(Tao tao,Vec X,PetscReal *f_reg,Vec G_reg,void *ptr) 88 { 89 AppCtx *user = (AppCtx*)ptr; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 /* compute regularizer objective 94 * f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x */ 95 ierr = VecCopy(X,user->workN2);CHKERRQ(ierr); 96 ierr = VecPow(user->workN2,2.);CHKERRQ(ierr); 97 ierr = VecShift(user->workN2,user->eps*user->eps);CHKERRQ(ierr); 98 ierr = VecSqrtAbs(user->workN2);CHKERRQ(ierr); 99 ierr = VecCopy(user->workN2, user->workN3);CHKERRQ(ierr); 100 ierr = VecShift(user->workN2,-user->eps);CHKERRQ(ierr); 101 ierr = VecSum(user->workN2,f_reg);CHKERRQ(ierr); 102 *f_reg *= user->lambda; 103 /* compute regularizer gradient = lambda*x */ 104 ierr = VecPointwiseDivide(G_reg,X,user->workN3);CHKERRQ(ierr); 105 ierr = VecScale(G_reg,user->lambda);CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108 109 /*------------------------------------------------------------*/ 110 111 PetscErrorCode RegularizerObjectiveAndGradient2(Tao tao,Vec X,PetscReal *f_reg,Vec G_reg,void *ptr) 112 { 113 AppCtx *user = (AppCtx*)ptr; 114 PetscErrorCode ierr; 115 PetscReal temp; 116 117 PetscFunctionBegin; 118 /* compute regularizer objective = lambda*|z|_2^2 */ 119 ierr = VecDot(X,X,&temp);CHKERRQ(ierr); 120 *f_reg = 0.5*user->lambda*temp; 121 /* compute regularizer gradient = lambda*z */ 122 ierr = VecCopy(X,G_reg);CHKERRQ(ierr); 123 ierr = VecScale(G_reg,user->lambda);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 /*------------------------------------------------------------*/ 128 129 static PetscErrorCode HessianMisfit(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr) 130 { 131 PetscFunctionBegin; 132 PetscFunctionReturn(0); 133 } 134 135 /*------------------------------------------------------------*/ 136 137 static PetscErrorCode HessianReg(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr) 138 { 139 AppCtx *user = (AppCtx*)ptr; 140 PetscErrorCode ierr; 141 142 PetscFunctionBegin; 143 ierr = MatMult(user->D,x,user->workN);CHKERRQ(ierr); 144 ierr = VecPow(user->workN2,2.);CHKERRQ(ierr); 145 ierr = VecShift(user->workN2,user->eps*user->eps);CHKERRQ(ierr); 146 ierr = VecSqrtAbs(user->workN2);CHKERRQ(ierr); 147 ierr = VecShift(user->workN2,-user->eps);CHKERRQ(ierr); 148 ierr = VecReciprocal(user->workN2);CHKERRQ(ierr); 149 ierr = VecScale(user->workN2,user->eps*user->eps);CHKERRQ(ierr); 150 ierr = MatDiagonalSet(H,user->workN2,INSERT_VALUES);CHKERRQ(ierr); 151 PetscFunctionReturn(0); 152 } 153 154 /*------------------------------------------------------------*/ 155 156 PetscErrorCode FullObjGrad(Tao tao,Vec X,PetscReal *f,Vec g,void *ptr) 157 { 158 AppCtx *user = (AppCtx*)ptr; 159 PetscErrorCode ierr; 160 PetscReal f_reg; 161 162 PetscFunctionBegin; 163 /* Objective 0.5*||Ax-b||_2^2 + lambda*||x||_2^2*/ 164 ierr = MatMult(user->A,X,user->workM);CHKERRQ(ierr); 165 ierr = VecAXPY(user->workM,-1,user->b);CHKERRQ(ierr); 166 ierr = VecDot(user->workM,user->workM,f);CHKERRQ(ierr); 167 ierr = VecNorm(X,NORM_2,&f_reg);CHKERRQ(ierr); 168 *f *= 0.5; 169 *f += user->lambda*f_reg*f_reg; 170 /* Gradient. ATAx-ATb + 2*lambda*x */ 171 ierr = MatMult(user->ATA,X,user->workN);CHKERRQ(ierr); 172 ierr = MatMultTranspose(user->A,user->b,user->workN2);CHKERRQ(ierr); 173 ierr = VecWAXPY(g,-1.,user->workN2,user->workN);CHKERRQ(ierr); 174 ierr = VecAXPY(g,2*user->lambda,X);CHKERRQ(ierr); 175 PetscFunctionReturn(0); 176 } 177 /*------------------------------------------------------------*/ 178 179 static PetscErrorCode HessianFull(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr) 180 { 181 PetscFunctionBegin; 182 PetscFunctionReturn(0); 183 } 184 /*------------------------------------------------------------*/ 185 186 PetscErrorCode InitializeUserData(AppCtx *user) 187 { 188 char dataFile[] = "tomographyData_A_b_xGT"; /* Matrix A and vectors b, xGT(ground truth) binary files generated by Matlab. Debug: change from "tomographyData_A_b_xGT" to "cs1Data_A_b_xGT". */ 189 PetscViewer fd; /* used to load data from file */ 190 PetscErrorCode ierr; 191 PetscInt k,n; 192 PetscScalar v; 193 PetscFunctionBegin; 194 195 /* Load the A matrix, b vector, and xGT vector from a binary file. */ 196 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,dataFile,FILE_MODE_READ,&fd);CHKERRQ(ierr); 197 ierr = MatCreate(PETSC_COMM_WORLD,&user->A);CHKERRQ(ierr); 198 ierr = MatSetType(user->A,MATAIJ);CHKERRQ(ierr); 199 ierr = MatLoad(user->A,fd);CHKERRQ(ierr); 200 ierr = VecCreate(PETSC_COMM_WORLD,&user->b);CHKERRQ(ierr); 201 ierr = VecLoad(user->b,fd);CHKERRQ(ierr); 202 ierr = VecCreate(PETSC_COMM_WORLD,&user->xGT);CHKERRQ(ierr); 203 ierr = VecLoad(user->xGT,fd);CHKERRQ(ierr); 204 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 205 206 ierr = MatGetSize(user->A,&user->M,&user->N);CHKERRQ(ierr); 207 208 ierr = MatCreate(PETSC_COMM_WORLD,&user->D);CHKERRQ(ierr); 209 ierr = MatSetSizes(user->D,PETSC_DECIDE,PETSC_DECIDE,user->N,user->N);CHKERRQ(ierr); 210 ierr = MatSetFromOptions(user->D);CHKERRQ(ierr); 211 ierr = MatSetUp(user->D);CHKERRQ(ierr); 212 for (k=0; k<user->N; k++) { 213 v = 1.0; 214 n = k+1; 215 if (k< user->N -1) { 216 ierr = MatSetValues(user->D,1,&k,1,&n,&v,INSERT_VALUES);CHKERRQ(ierr); 217 } 218 v = -1.0; 219 ierr = MatSetValues(user->D,1,&k,1,&k,&v,INSERT_VALUES);CHKERRQ(ierr); 220 } 221 ierr = MatAssemblyBegin(user->D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 222 ierr = MatAssemblyEnd(user->D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 223 224 ierr = MatTransposeMatMult(user->D,user->D,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&user->DTD);CHKERRQ(ierr); 225 226 ierr = MatCreate(PETSC_COMM_WORLD,&user->Hz);CHKERRQ(ierr); 227 ierr = MatSetSizes(user->Hz,PETSC_DECIDE,PETSC_DECIDE,user->N,user->N);CHKERRQ(ierr); 228 ierr = MatSetFromOptions(user->Hz);CHKERRQ(ierr); 229 ierr = MatSetUp(user->Hz);CHKERRQ(ierr); 230 ierr = MatAssemblyBegin(user->Hz,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 231 ierr = MatAssemblyEnd(user->Hz,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 232 233 ierr = VecCreate(PETSC_COMM_WORLD,&(user->x));CHKERRQ(ierr); 234 ierr = VecCreate(PETSC_COMM_WORLD,&(user->workM));CHKERRQ(ierr); 235 ierr = VecCreate(PETSC_COMM_WORLD,&(user->workN));CHKERRQ(ierr); 236 ierr = VecCreate(PETSC_COMM_WORLD,&(user->workN2));CHKERRQ(ierr); 237 ierr = VecSetSizes(user->x,PETSC_DECIDE,user->N);CHKERRQ(ierr); 238 ierr = VecSetSizes(user->workM,PETSC_DECIDE,user->M);CHKERRQ(ierr); 239 ierr = VecSetSizes(user->workN,PETSC_DECIDE,user->N);CHKERRQ(ierr); 240 ierr = VecSetSizes(user->workN2,PETSC_DECIDE,user->N);CHKERRQ(ierr); 241 ierr = VecSetFromOptions(user->x);CHKERRQ(ierr); 242 ierr = VecSetFromOptions(user->workM);CHKERRQ(ierr); 243 ierr = VecSetFromOptions(user->workN);CHKERRQ(ierr); 244 ierr = VecSetFromOptions(user->workN2);CHKERRQ(ierr); 245 246 ierr = VecDuplicate(user->workN,&(user->workN3));CHKERRQ(ierr); 247 ierr = VecDuplicate(user->x,&(user->xlb));CHKERRQ(ierr); 248 ierr = VecDuplicate(user->x,&(user->xub));CHKERRQ(ierr); 249 ierr = VecDuplicate(user->x,&(user->c));CHKERRQ(ierr); 250 ierr = VecSet(user->xlb,0.0);CHKERRQ(ierr); 251 ierr = VecSet(user->c,0.0);CHKERRQ(ierr); 252 ierr = VecSet(user->xub,PETSC_INFINITY);CHKERRQ(ierr); 253 254 ierr = MatTransposeMatMult(user->A,user->A, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &(user->ATA));CHKERRQ(ierr); 255 ierr = MatTransposeMatMult(user->A,user->A, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &(user->Hx));CHKERRQ(ierr); 256 ierr = MatTransposeMatMult(user->A,user->A, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &(user->HF));CHKERRQ(ierr); 257 258 ierr = MatAssemblyBegin(user->ATA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 259 ierr = MatAssemblyEnd(user->ATA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 260 ierr = MatAssemblyBegin(user->Hx,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 261 ierr = MatAssemblyEnd(user->Hx,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 262 ierr = MatAssemblyBegin(user->HF,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 263 ierr = MatAssemblyEnd(user->HF,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 264 265 user->lambda = 1.e-8; 266 user->eps = 1.e-3; 267 user->reg = 2; 268 user->mumin = 5.e-6; 269 270 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Configure separable objection example", "tomographyADMM.c");CHKERRQ(ierr); 271 ierr = PetscOptionsInt("-reg","Regularization scheme for z solver (1,2)", "tomographyADMM.c", user->reg, &(user->reg), NULL);CHKERRQ(ierr); 272 ierr = PetscOptionsReal("-lambda", "The regularization multiplier. 1 default", "tomographyADMM.c", user->lambda, &(user->lambda), NULL);CHKERRQ(ierr); 273 ierr = PetscOptionsReal("-eps", "L1 norm epsilon padding", "tomographyADMM.c", user->eps, &(user->eps), NULL);CHKERRQ(ierr); 274 ierr = PetscOptionsReal("-mumin", "Minimum value for ADMM spectral penalty", "tomographyADMM.c", user->mumin, &(user->mumin), NULL);CHKERRQ(ierr); 275 ierr = PetscOptionsEnd();CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278 279 /*------------------------------------------------------------*/ 280 281 PetscErrorCode DestroyContext(AppCtx *user) 282 { 283 PetscErrorCode ierr; 284 285 PetscFunctionBegin; 286 ierr = MatDestroy(&user->A);CHKERRQ(ierr); 287 ierr = MatDestroy(&user->ATA);CHKERRQ(ierr); 288 ierr = MatDestroy(&user->Hx);CHKERRQ(ierr); 289 ierr = MatDestroy(&user->Hz);CHKERRQ(ierr); 290 ierr = MatDestroy(&user->HF);CHKERRQ(ierr); 291 ierr = MatDestroy(&user->D);CHKERRQ(ierr); 292 ierr = MatDestroy(&user->DTD);CHKERRQ(ierr); 293 ierr = VecDestroy(&user->xGT);CHKERRQ(ierr); 294 ierr = VecDestroy(&user->xlb);CHKERRQ(ierr); 295 ierr = VecDestroy(&user->xub);CHKERRQ(ierr); 296 ierr = VecDestroy(&user->b);CHKERRQ(ierr); 297 ierr = VecDestroy(&user->x);CHKERRQ(ierr); 298 ierr = VecDestroy(&user->c);CHKERRQ(ierr); 299 ierr = VecDestroy(&user->workN3);CHKERRQ(ierr); 300 ierr = VecDestroy(&user->workN2);CHKERRQ(ierr); 301 ierr = VecDestroy(&user->workN);CHKERRQ(ierr); 302 ierr = VecDestroy(&user->workM);CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 /*------------------------------------------------------------*/ 307 308 int main(int argc,char **argv) 309 { 310 PetscErrorCode ierr; 311 Tao tao,misfit,reg; 312 PetscReal v1,v2; 313 AppCtx* user; 314 PetscViewer fd; 315 char resultFile[] = "tomographyResult_x"; 316 317 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 318 ierr = PetscNew(&user);CHKERRQ(ierr); 319 ierr = InitializeUserData(user);CHKERRQ(ierr); 320 321 ierr = TaoCreate(PETSC_COMM_WORLD, &tao);CHKERRQ(ierr); 322 ierr = TaoSetType(tao, TAOADMM);CHKERRQ(ierr); 323 ierr = TaoSetInitialVector(tao, user->x);CHKERRQ(ierr); 324 /* f(x) + g(x) for parent tao */ 325 ierr = TaoADMMSetSpectralPenalty(tao,1.);CHKERRQ(ierr); 326 ierr = TaoSetObjectiveAndGradientRoutine(tao, FullObjGrad, (void*)user);CHKERRQ(ierr); 327 ierr = MatShift(user->HF,user->lambda);CHKERRQ(ierr); 328 ierr = TaoSetHessianRoutine(tao, user->HF, user->HF, HessianFull, (void*)user);CHKERRQ(ierr); 329 330 /* f(x) for misfit tao */ 331 ierr = TaoADMMSetMisfitObjectiveAndGradientRoutine(tao, MisfitObjectiveAndGradient, (void*)user);CHKERRQ(ierr); 332 ierr = TaoADMMSetMisfitHessianRoutine(tao, user->Hx, user->Hx, HessianMisfit, (void*)user);CHKERRQ(ierr); 333 ierr = TaoADMMSetMisfitHessianChangeStatus(tao,PETSC_FALSE);CHKERRQ(ierr); 334 ierr = TaoADMMSetMisfitConstraintJacobian(tao,user->D,user->D,NullJacobian,(void*)user);CHKERRQ(ierr); 335 336 /* g(x) for regularizer tao */ 337 if (user->reg == 1) { 338 ierr = TaoADMMSetRegularizerObjectiveAndGradientRoutine(tao, RegularizerObjectiveAndGradient1, (void*)user);CHKERRQ(ierr); 339 ierr = TaoADMMSetRegularizerHessianRoutine(tao, user->Hz, user->Hz, HessianReg, (void*)user);CHKERRQ(ierr); 340 ierr = TaoADMMSetRegHessianChangeStatus(tao,PETSC_TRUE);CHKERRQ(ierr); 341 } else if (user->reg == 2) { 342 ierr = TaoADMMSetRegularizerObjectiveAndGradientRoutine(tao, RegularizerObjectiveAndGradient2, (void*)user);CHKERRQ(ierr); 343 ierr = MatShift(user->Hz,1);CHKERRQ(ierr); 344 ierr = MatScale(user->Hz,user->lambda);CHKERRQ(ierr); 345 ierr = TaoADMMSetRegularizerHessianRoutine(tao, user->Hz, user->Hz, HessianMisfit, (void*)user);CHKERRQ(ierr); 346 ierr = TaoADMMSetRegHessianChangeStatus(tao,PETSC_TRUE);CHKERRQ(ierr); 347 } else if (user->reg != 3) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_UNKNOWN_TYPE, "Incorrect Reg type"); /* TaoShell case */ 348 349 /* Set type for the misfit solver */ 350 ierr = TaoADMMGetMisfitSubsolver(tao, &misfit);CHKERRQ(ierr); 351 ierr = TaoADMMGetRegularizationSubsolver(tao, ®);CHKERRQ(ierr); 352 ierr = TaoSetType(misfit,TAONLS);CHKERRQ(ierr); 353 if (user->reg == 3) { 354 ierr = TaoSetType(reg,TAOSHELL);CHKERRQ(ierr); 355 ierr = TaoShellSetContext(reg, (void*) user);CHKERRQ(ierr); 356 ierr = TaoShellSetSolve(reg, TaoShellSolve_SoftThreshold);CHKERRQ(ierr); 357 } else { 358 ierr = TaoSetType(reg,TAONLS);CHKERRQ(ierr); 359 } 360 ierr = TaoSetVariableBounds(misfit,user->xlb,user->xub);CHKERRQ(ierr); 361 362 /* Soft Thresholding solves the ADMM problem with the L1 regularizer lambda*||z||_1 and the x-z=0 constraint */ 363 ierr = TaoADMMSetRegularizerCoefficient(tao, user->lambda);CHKERRQ(ierr); 364 ierr = TaoADMMSetRegularizerConstraintJacobian(tao,NULL,NULL,NullJacobian,(void*)user);CHKERRQ(ierr); 365 ierr = TaoADMMSetMinimumSpectralPenalty(tao,user->mumin);CHKERRQ(ierr); 366 367 ierr = TaoADMMSetConstraintVectorRHS(tao,user->c);CHKERRQ(ierr); 368 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 369 ierr = TaoSolve(tao);CHKERRQ(ierr); 370 371 /* Save x (reconstruction of object) vector to a binary file, which maybe read from Matlab and convert to a 2D image for comparison. */ 372 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,resultFile,FILE_MODE_WRITE,&fd);CHKERRQ(ierr); 373 ierr = VecView(user->x,fd);CHKERRQ(ierr); 374 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 375 376 /* compute the error */ 377 ierr = VecAXPY(user->x,-1,user->xGT);CHKERRQ(ierr); 378 ierr = VecNorm(user->x,NORM_2,&v1);CHKERRQ(ierr); 379 ierr = VecNorm(user->xGT,NORM_2,&v2);CHKERRQ(ierr); 380 ierr = PetscPrintf(PETSC_COMM_WORLD, "relative reconstruction error: ||x-xGT||/||xGT|| = %6.4e.\n", (double)(v1/v2));CHKERRQ(ierr); 381 382 /* Free TAO data structures */ 383 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 384 ierr = DestroyContext(user);CHKERRQ(ierr); 385 ierr = PetscFree(user);CHKERRQ(ierr); 386 ierr = PetscFinalize(); 387 return ierr; 388 } 389 390 /*TEST 391 392 build: 393 requires: !complex !single !__float128 !defined(PETSC_USE_64BIT_INDICES) 394 395 test: 396 suffix: 1 397 localrunfiles: tomographyData_A_b_xGT 398 args: -lambda 1.e-8 -tao_monitor -tao_type nls -tao_nls_pc_type icc 399 400 test: 401 suffix: 2 402 localrunfiles: tomographyData_A_b_xGT 403 args: -reg 2 -lambda 1.e-8 -tao_admm_dual_update update_basic -tao_admm_regularizer_type regularizer_user -tao_max_it 20 -tao_monitor -tao_admm_tolerance_update_factor 1.e-8 -misfit_tao_nls_pc_type icc -misfit_tao_monitor -reg_tao_monitor 404 405 test: 406 suffix: 3 407 localrunfiles: tomographyData_A_b_xGT 408 args: -lambda 1.e-8 -tao_admm_dual_update update_basic -tao_admm_regularizer_type regularizer_soft_thresh -tao_max_it 20 -tao_monitor -tao_admm_tolerance_update_factor 1.e-8 -misfit_tao_nls_pc_type icc -misfit_tao_monitor 409 410 test: 411 suffix: 4 412 localrunfiles: tomographyData_A_b_xGT 413 args: -lambda 1.e-8 -tao_admm_dual_update update_adaptive -tao_admm_regularizer_type regularizer_soft_thresh -tao_max_it 20 -tao_monitor -misfit_tao_monitor -misfit_tao_nls_pc_type icc 414 415 test: 416 suffix: 5 417 localrunfiles: tomographyData_A_b_xGT 418 args: -reg 2 -lambda 1.e-8 -tao_admm_dual_update update_adaptive -tao_admm_regularizer_type regularizer_user -tao_max_it 20 -tao_monitor -tao_admm_tolerance_update_factor 1.e-8 -misfit_tao_monitor -reg_tao_monitor -misfit_tao_nls_pc_type icc 419 420 test: 421 suffix: 6 422 localrunfiles: tomographyData_A_b_xGT 423 args: -reg 3 -lambda 1.e-8 -tao_admm_dual_update update_adaptive -tao_admm_regularizer_type regularizer_user -tao_max_it 20 -tao_monitor -tao_admm_tolerance_update_factor 1.e-8 -misfit_tao_monitor -reg_tao_monitor -misfit_tao_nls_pc_type icc 424 425 TEST*/ 426