1 #include "petsc.h" 2 #include "petscblaslapack.h" 3 #include "taolapack.h" 4 5 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "estsv" 9 static PetscErrorCode estsv(PetscInt n, PetscReal *r, PetscInt ldr, PetscReal *svmin, PetscReal *z) { 10 11 PetscBLASInt blas1=1, blasn=n, blasnmi, blasj, blasldr = ldr; 12 PetscInt i,j; 13 PetscReal e,temp,w,wm,ynorm,znorm,s,sm; 14 PetscFunctionBegin; 15 for (i=0;i<n;i++) { 16 z[i]=0.0; 17 } 18 19 e = PetscAbs(r[0]); 20 if (e == 0.0) { 21 *svmin = 0.0; 22 z[0] = 1.0; 23 } else { 24 /* Solve R'*y = e */ 25 for (i=0;i<n;i++) { 26 27 /* Scale y. The scaling factor (0.01) reduces the number of scalings */ 28 if (z[i] >= 0.0) 29 e=-PetscAbs(e); 30 else 31 e = PetscAbs(e); 32 33 if (PetscAbs(e - z[i]) > PetscAbs(r[i + ldr*i])) { 34 temp = PetscMin(0.01,PetscAbs(r[i + ldr*i]))/ 35 PetscAbs(e-z[i]); 36 BLASscal_(&blasn, &temp, z, &blas1); 37 e = temp*e; 38 } 39 40 /* Determine the two possible choices of y[i] */ 41 if (r[i + ldr*i] == 0.0) 42 w = wm = 1.0; 43 else { 44 w = (e - z[i]) / r[i + ldr*i]; 45 wm = - (e + z[i]) / r[i + ldr*i]; 46 } 47 48 /* Chose y[i] based on the predicted value of y[j] for j>i */ 49 s = PetscAbs(e - z[i]); 50 sm = PetscAbs(e + z[i]); 51 for (j=i+1;j<n;j++) { 52 sm += PetscAbs(z[j] + wm * r[i + ldr*j]); 53 } 54 if (i < n-1) { 55 blasnmi = n-i-1; 56 BLASaxpy_(&blasnmi, &w, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1); 57 s += BLASasum_(&blasnmi, &z[i+1], &blas1); 58 } 59 if (s < sm) { 60 temp = wm - w; 61 w = wm; 62 if (i < n-1) { 63 BLASaxpy_(&blasnmi, &temp, &r[i + ldr*(i+1)], &blasldr, &z[i+1], &blas1); 64 } 65 } 66 z[i] = w; 67 } 68 69 ynorm = BLASnrm2_(&blasn, z, &blas1); 70 71 /* Solve R*z = y */ 72 for (j=n-1; j>=0; j--) { 73 74 /* Scale z */ 75 76 if (PetscAbs(z[j]) > PetscAbs(r[j + ldr*j])) { 77 temp = PetscMin(0.01, PetscAbs(r[j + ldr*j] / z[j])); 78 BLASscal_(&blasn, &temp, z, &blas1); 79 ynorm *=temp; 80 } 81 if (r[j + ldr*j] == 0) { 82 z[j] = 1.0; 83 } else { 84 z[j] = z[j] / r[j + ldr*j]; 85 } 86 temp = -z[j]; 87 blasj=j; 88 BLASaxpy_(&blasj,&temp,&r[0+ldr*j],&blas1,z,&blas1); 89 } 90 91 /* Compute svmin and normalize z */ 92 znorm = 1.0 / BLASnrm2_(&blasn, z, &blas1); 93 *svmin = ynorm*znorm; 94 BLASscal_(&blasn, &znorm, z, &blas1); 95 96 } 97 98 PetscFunctionReturn(0); 99 } 100 101 102 103 /* 104 c *********** 105 c 106 c Subroutine dgqt 107 c 108 c Given an n by n symmetric matrix A, an n-vector b, and a 109 c positive number delta, this subroutine determines a vector 110 c x which approximately minimizes the quadratic function 111 c 112 c f(x) = (1/2)*x'*A*x + b'*x 113 c 114 c subject to the Euclidean norm constraint 115 c 116 c norm(x) <= delta. 117 c 118 c This subroutine computes an approximation x and a Lagrange 119 c multiplier par such that either par is zero and 120 c 121 c norm(x) <= (1+rtol)*delta, 122 c 123 c or par is positive and 124 c 125 c abs(norm(x) - delta) <= rtol*delta. 126 c 127 c If xsol is the solution to the problem, the approximation x 128 c satisfies 129 c 130 c f(x) <= ((1 - rtol)**2)*f(xsol) 131 c 132 c The subroutine statement is 133 c 134 c subroutine dgqt(n,a,lda,b,delta,rtol,atol,itmax, 135 c par,f,x,info,z,wa1,wa2) 136 c 137 c where 138 c 139 c n is an integer variable. 140 c On entry n is the order of A. 141 c On exit n is unchanged. 142 c 143 c a is a double precision array of dimension (lda,n). 144 c On entry the full upper triangle of a must contain the 145 c full upper triangle of the symmetric matrix A. 146 c On exit the array contains the matrix A. 147 c 148 c lda is an integer variable. 149 c On entry lda is the leading dimension of the array a. 150 c On exit lda is unchanged. 151 c 152 c b is an double precision array of dimension n. 153 c On entry b specifies the linear term in the quadratic. 154 c On exit b is unchanged. 155 c 156 c delta is a double precision variable. 157 c On entry delta is a bound on the Euclidean norm of x. 158 c On exit delta is unchanged. 159 c 160 c rtol is a double precision variable. 161 c On entry rtol is the relative accuracy desired in the 162 c solution. Convergence occurs if 163 c 164 c f(x) <= ((1 - rtol)**2)*f(xsol) 165 c 166 c On exit rtol is unchanged. 167 c 168 c atol is a double precision variable. 169 c On entry atol is the absolute accuracy desired in the 170 c solution. Convergence occurs when 171 c 172 c norm(x) <= (1 + rtol)*delta 173 c 174 c max(-f(x),-f(xsol)) <= atol 175 c 176 c On exit atol is unchanged. 177 c 178 c itmax is an integer variable. 179 c On entry itmax specifies the maximum number of iterations. 180 c On exit itmax is unchanged. 181 c 182 c par is a double precision variable. 183 c On entry par is an initial estimate of the Lagrange 184 c multiplier for the constraint norm(x) <= delta. 185 c On exit par contains the final estimate of the multiplier. 186 c 187 c f is a double precision variable. 188 c On entry f need not be specified. 189 c On exit f is set to f(x) at the output x. 190 c 191 c x is a double precision array of dimension n. 192 c On entry x need not be specified. 193 c On exit x is set to the final estimate of the solution. 194 c 195 c info is an integer variable. 196 c On entry info need not be specified. 197 c On exit info is set as follows: 198 c 199 c info = 1 The function value f(x) has the relative 200 c accuracy specified by rtol. 201 c 202 c info = 2 The function value f(x) has the absolute 203 c accuracy specified by atol. 204 c 205 c info = 3 Rounding errors prevent further progress. 206 c On exit x is the best available approximation. 207 c 208 c info = 4 Failure to converge after itmax iterations. 209 c On exit x is the best available approximation. 210 c 211 c z is a double precision work array of dimension n. 212 c 213 c wa1 is a double precision work array of dimension n. 214 c 215 c wa2 is a double precision work array of dimension n. 216 c 217 c Subprograms called 218 c 219 c MINPACK-2 ...... destsv 220 c 221 c LAPACK ......... dpotrf 222 c 223 c Level 1 BLAS ... daxpy, dcopy, ddot, dnrm2, dscal 224 c 225 c Level 2 BLAS ... dtrmv, dtrsv 226 c 227 c MINPACK-2 Project. October 1993. 228 c Argonne National Laboratory and University of Minnesota. 229 c Brett M. Averick, Richard Carter, and Jorge J. More' 230 c 231 c *********** 232 */ 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "gqt" 236 PetscErrorCode gqt(PetscInt n, PetscReal *a, PetscInt lda, PetscReal *b, 237 PetscReal delta, PetscReal rtol, PetscReal atol, 238 PetscInt itmax, PetscReal *retpar, PetscReal *retf, 239 PetscReal *x, PetscInt *retinfo, PetscInt *retits, 240 PetscReal *z, PetscReal *wa1, PetscReal *wa2) 241 { 242 PetscErrorCode ierr; 243 PetscReal f=0.0,p001=0.001,p5=0.5,minusone=-1,delta2=delta*delta; 244 PetscInt iter, j, rednc,info; 245 PetscBLASInt indef; 246 PetscBLASInt blas1=1, blasn=n, iblas, blaslda = lda,blasldap1=lda+1,blasinfo; 247 PetscReal alpha, anorm, bnorm, parc, parf, parl, pars, par=*retpar, 248 paru, prod, rxnorm, rznorm=0.0, temp, xnorm; 249 250 PetscFunctionBegin; 251 iter = 0; 252 parf = 0.0; 253 xnorm = 0.0; 254 rxnorm = 0.0; 255 rednc = 0; 256 for (j=0; j<n; j++) { 257 x[j] = 0.0; 258 z[j] = 0.0; 259 } 260 261 /* Copy the diagonal and save A in its lower triangle */ 262 BLAScopy_(&blasn,a,&blasldap1, wa1, &blas1); 263 CHKMEMQ; 264 for (j=0;j<n-1;j++) { 265 iblas = n - j - 1; 266 BLAScopy_(&iblas,&a[j + lda*(j+1)], &blaslda, &a[j+1 + lda*j], &blas1); 267 CHKMEMQ; 268 } 269 270 /* Calculate the l1-norm of A, the Gershgorin row sums, and the 271 l2-norm of b */ 272 anorm = 0.0; 273 for (j=0;j<n;j++) { 274 wa2[j] = BLASasum_(&blasn, &a[0 + lda*j], &blas1); 275 CHKMEMQ; 276 anorm = PetscMax(anorm,wa2[j]); 277 } 278 for (j=0;j<n;j++) { 279 wa2[j] = wa2[j] - PetscAbs(wa1[j]); 280 } 281 bnorm = BLASnrm2_(&blasn,b,&blas1); 282 CHKMEMQ; 283 /* Calculate a lower bound, pars, for the domain of the problem. 284 Also calculate an upper bound, paru, and a lower bound, parl, 285 for the Lagrange multiplier. */ 286 pars = parl = paru = -anorm; 287 for (j=0;j<n;j++) { 288 pars = PetscMax(pars, -wa1[j]); 289 parl = PetscMax(parl, wa1[j] + wa2[j]); 290 paru = PetscMax(paru, -wa1[j] + wa2[j]); 291 } 292 parl = PetscMax(bnorm/delta - parl,pars); 293 parl = PetscMax(0.0,parl); 294 paru = PetscMax(0.0, bnorm/delta + paru); 295 296 /* If the input par lies outside of the interval (parl, paru), 297 set par to the closer endpoint. */ 298 299 par = PetscMax(par,parl); 300 par = PetscMin(par,paru); 301 302 /* Special case: parl == paru */ 303 paru = PetscMax(paru, (1.0 + rtol)*parl); 304 305 /* Beginning of an iteration */ 306 307 info = 0; 308 for (iter=1;iter<=itmax;iter++) { 309 310 /* Safeguard par */ 311 if (par <= pars && paru > 0) { 312 par = PetscMax(p001, PetscSqrtScalar(parl/paru)) * paru; 313 } 314 315 /* Copy the lower triangle of A into its upper triangle and 316 compute A + par*I */ 317 318 for (j=0;j<n-1;j++) { 319 iblas = n - j - 1; 320 BLAScopy_(&iblas,&a[j+1 + j*lda], &blas1, 321 &a[j + (j+1)*lda], &blaslda); 322 CHKMEMQ; 323 } 324 for (j=0;j<n;j++) { 325 a[j + j*lda] = wa1[j] + par; 326 } 327 328 /* Attempt the Cholesky factorization of A without referencing 329 the lower triangular part. */ 330 LAPACKpotrf_("U",&blasn,a,&blaslda,&indef); 331 CHKMEMQ; 332 333 /* Case 1: A + par*I is pos. def. */ 334 if (indef == 0) { 335 336 /* Compute an approximate solution x and save the 337 last value of par with A + par*I pos. def. */ 338 339 parf = par; 340 BLAScopy_(&blasn, b, &blas1, wa2, &blas1); 341 CHKMEMQ; 342 LAPACKtrtrs_("U","T","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo); 343 CHKMEMQ; 344 rxnorm = BLASnrm2_(&blasn, wa2, &blas1); 345 LAPACKtrtrs_("U","N","N",&blasn,&blas1,a,&blaslda,wa2,&blasn,&blasinfo); 346 CHKMEMQ; 347 BLAScopy_(&blasn, wa2, &blas1, x, &blas1); 348 CHKMEMQ; 349 BLASscal_(&blasn, &minusone, x, &blas1); 350 CHKMEMQ; 351 xnorm = BLASnrm2_(&blasn, x, &blas1); 352 CHKMEMQ; 353 354 /* Test for convergence */ 355 if (PetscAbs(xnorm - delta) <= rtol*delta || 356 (par == 0 && xnorm <= (1.0+rtol)*delta)) { 357 info = 1; 358 } 359 360 /* Compute a direction of negative curvature and use this 361 information to improve pars. */ 362 363 iblas=blasn*blasn; 364 365 ierr = estsv(n,a,lda,&rznorm,z); CHKERRQ(ierr); 366 CHKMEMQ; 367 pars = PetscMax(pars, par-rznorm*rznorm); 368 369 /* Compute a negative curvature solution of the form 370 x + alpha*z, where norm(x+alpha*z)==delta */ 371 372 rednc = 0; 373 if (xnorm < delta) { 374 375 /* Compute alpha */ 376 prod = BLASdot_(&blasn, z, &blas1, x, &blas1) / delta; 377 temp = (delta - xnorm)*((delta + xnorm)/delta); 378 alpha = temp/(PetscAbs(prod) + PetscSqrtScalar(prod*prod + temp/delta)); 379 if (prod >= 0) 380 alpha = PetscAbs(alpha); 381 else 382 alpha =-PetscAbs(alpha); 383 384 385 /* Test to decide if the negative curvature step 386 produces a larger reduction than with z=0 */ 387 388 rznorm = PetscAbs(alpha) * rznorm; 389 if ((rznorm*rznorm + par*xnorm*xnorm)/(delta2) <= par) { 390 rednc = 1; 391 } 392 393 /* Test for convergence */ 394 if (p5 * rznorm*rznorm / delta2 <= 395 rtol*(1.0-p5*rtol)*(par + rxnorm*rxnorm/delta2)) { 396 info = 1; 397 } else if (info == 0 && 398 (p5*(par + rxnorm*rxnorm/delta2) <= atol/delta2)) { 399 info = 2; 400 } 401 } 402 403 /* Compute the Newton correction parc to par. */ 404 if (xnorm == 0) { 405 parc = -par; 406 } else { 407 BLAScopy_(&blasn, x, &blas1, wa2, &blas1); 408 CHKMEMQ; 409 temp = 1.0/xnorm; 410 BLASscal_(&blasn, &temp, wa2, &blas1); 411 CHKMEMQ; 412 LAPACKtrtrs_("U","T","N",&blasn, &blas1, a, &blaslda, wa2, &blasn, &blasinfo); 413 CHKMEMQ; 414 temp = BLASnrm2_(&blasn, wa2, &blas1); 415 parc = (xnorm - delta)/(delta*temp*temp); 416 } 417 418 /* update parl or paru */ 419 if (xnorm > delta) { 420 parl = PetscMax(parl, par); 421 } else if (xnorm < delta) { 422 paru = PetscMin(paru, par); 423 } 424 } else { 425 /* Case 2: A + par*I is not pos. def. */ 426 427 /* Use the rank information from the Cholesky 428 decomposition to update par. */ 429 430 if (indef > 1) { 431 432 /* Restore column indef to A + par*I. */ 433 iblas = indef - 1; 434 BLAScopy_(&iblas,&a[indef-1 + 0*lda],&blaslda, 435 &a[0 + (indef-1)*lda],&blas1); 436 CHKMEMQ; 437 a[indef-1 + (indef-1)*lda] = wa1[indef-1] + par; 438 439 /* compute parc. */ 440 441 BLAScopy_(&iblas,&a[0 + (indef-1)*lda], &blas1, wa2, &blas1); 442 CHKMEMQ; 443 LAPACKtrtrs_("U","T","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo); 444 CHKMEMQ; 445 BLAScopy_(&iblas,wa2,&blas1,&a[0 + (indef-1)*lda],&blas1); 446 CHKMEMQ; 447 temp = BLASnrm2_(&iblas,&a[0 + (indef-1)*lda],&blas1); 448 CHKMEMQ; 449 a[indef-1 + (indef-1)*lda] -= temp*temp; 450 LAPACKtrtrs_("U","N","N",&iblas,&blas1,a,&blaslda,wa2,&blasn,&blasinfo); 451 CHKMEMQ; 452 } 453 454 wa2[indef-1] = -1.0; 455 iblas = indef; 456 temp = BLASnrm2_(&iblas,wa2,&blas1); 457 parc = - a[indef-1 + (indef-1)*lda]/(temp*temp); 458 pars = PetscMax(pars,par+parc); 459 460 /* If necessary, increase paru slightly. 461 This is needed because in some exceptional situations 462 paru is the optimal value of par. */ 463 464 paru = PetscMax(paru, (1.0+rtol)*pars); 465 } 466 467 /* Use pars to update parl */ 468 parl = PetscMax(parl,pars); 469 470 /* Test for termination. */ 471 if (info == 0) { 472 if (iter == itmax) info=4; 473 if (paru <= (1.0+p5*rtol)*pars) info=3; 474 if (paru == 0.0) info = 2; 475 } 476 477 /* If exiting, store the best approximation and restore 478 the upper triangle of A. */ 479 480 if (info != 0) { 481 482 /* Compute the best current estimates for x and f. */ 483 484 par = parf; 485 f = -p5 * (rxnorm*rxnorm + par*xnorm*xnorm); 486 if (rednc) { 487 f = -p5 * (rxnorm*rxnorm + par*delta*delta - rznorm*rznorm); 488 BLASaxpy_(&blasn, &alpha, z, &blas1, x, &blas1); 489 CHKMEMQ; 490 } 491 492 /* Restore the upper triangle of A */ 493 494 for (j = 0; j<n; j++) { 495 iblas = n - j - 1; 496 BLAScopy_(&iblas,&a[j+1 + j*lda],&blas1, &a[j + (j+1)*lda],&blaslda); 497 CHKMEMQ; 498 } 499 iblas = lda+1; 500 BLAScopy_(&blasn,wa1,&blas1,a,&iblas); 501 CHKMEMQ; 502 break; 503 } 504 par = PetscMax(parl,par+parc); 505 } 506 *retpar = par; 507 *retf = f; 508 *retinfo = info; 509 *retits = iter; 510 CHKMEMQ; 511 PetscFunctionReturn(0); 512 } 513