Actual source code: lmvmmat.c
1: #include "lmvmmat.h"
2: #include "tao_general.h"
3: #include "tao_solver.h"
4: #include "taovec.h"
6: #define Scale_None 0
7: #define Scale_Scalar 1
8: #define Scale_Broyden 2
9: #define Scale_Types 3
11: #define Rescale_None 0
12: #define Rescale_Scalar 1
13: #define Rescale_GL 2
14: #define Rescale_Types 3
16: #define Limit_None 0
17: #define Limit_Average 1
18: #define Limit_Relative 2
19: #define Limit_Absolute 3
20: #define Limit_Types 4
22: #define TAO_ZER_SAFEGUARD 1e-8
23: #define TAO_INF_SAFEGUARD 1e+8
25: static const char *Scale_Table[64] = {
26: "none", "scalar", "broyden"
27: };
29: static const char *Rescale_Table[64] = {
30: "none", "scalar", "gl"
31: };
33: static const char *Limit_Table[64] = {
34: "none", "average", "relative", "absolute"
35: };
39: TaoLMVMMat::TaoLMVMMat(TaoVec *tv)
40: {
41: // Set default values
42: lm = 5;
43: eps = 0.0;
45: limitType = Limit_None;
46: scaleType = Scale_Broyden;
47: rScaleType = Rescale_Scalar;
49: s_alpha = 1.0;
50: r_alpha = 1.0;
51: r_beta = 0.5;
52: mu = 1.0;
53: nu = 100.0;
55: phi = 0.125;
57: scalar_history = 1;
58: rescale_history = 1;
60: delta_min = 1e-7;
61: delta_max = 100;
63: // Begin configuration
64: TaoOptionsHead("Limited memory matrix options");
65: TaoOptionInt("-tao_lmm_vectors", "vectors to use for approximation", "", lm, &lm, 0);
66: TaoOptionDouble("-tao_lmm_limit_mu", "mu limiting factor", "", mu, &mu, 0);
67: TaoOptionDouble("-tao_lmm_limit_nu", "nu limiting factor", "", nu, &nu, 0);
68: TaoOptionDouble("-tao_lmm_broyden_phi", "phi factor for Broyden scaling", "", phi, &phi, 0);
69: TaoOptionDouble("-tao_lmm_scalar_alpha", "alpha factor for scalar scaling", "", s_alpha, &s_alpha, 0);
70: TaoOptionDouble("-tao_lmm_rescale_alpha", "alpha factor for rescaling diagonal", "", r_alpha, &r_alpha, 0);
71: TaoOptionDouble("-tao_lmm_rescale_beta", "beta factor for rescaling diagonal", "", r_beta, &r_beta, 0);
72: TaoOptionInt("-tao_lmm_scalar_history", "amount of history for scalar scaling", "", scalar_history, &scalar_history, 0);
73: TaoOptionInt("-tao_lmm_rescale_history", "amount of history for rescaling diagonal", "", rescale_history, &rescale_history, 0);
74: TaoOptionDouble("-tao_lmm_eps", "rejection tolerance", "", eps, &eps, 0);
75: TaoOptionList("-tao_lmm_scale_type", "scale type", "", Scale_Table, Scale_Types, Scale_Table[scaleType], &scaleType, 0);
76: TaoOptionList("-tao_lmm_rescale_type", "rescale type", "", Rescale_Table, Rescale_Types, Rescale_Table[rScaleType], &rScaleType, 0);
77: TaoOptionList("-tao_lmm_limit_type", "limit type", "", Limit_Table, Limit_Types, Limit_Table[limitType], &limitType, 0);
78: TaoOptionDouble("-tao_lmm_delta_min", "minimum delta value", "", delta_min, &delta_min, 0);
79: TaoOptionDouble("-tao_lmm_delta_max", "maximum delta value", "", delta_max, &delta_max, 0);
80: TaoOptionsTail();
82: // Complete configuration
83: rescale_history = TaoMin(rescale_history, lm);
85: // Perform allocations
86: tv->CloneVecs(lm+1, &S);
87: tv->CloneVecs(lm+1, &Y);
88: tv->Clone(&D);
89: tv->Clone(&U);
90: tv->Clone(&V);
91: tv->Clone(&W);
92: tv->Clone(&P);
93: tv->Clone(&Q);
95: rho = new double[lm+1];
96: beta = new double[lm+1];
98: yy_history = new double[TaoMax(scalar_history, 1)];
99: ys_history = new double[TaoMax(scalar_history, 1)];
100: ss_history = new double[TaoMax(scalar_history, 1)];
102: yy_rhistory = new double[TaoMax(rescale_history, 1)];
103: ys_rhistory = new double[TaoMax(rescale_history, 1)];
104: ss_rhistory = new double[TaoMax(rescale_history, 1)];
106: // Finish initializations
107: lmnow = 0;
108: iter = 0;
109: updates = 0;
110: rejects = 0;
111: delta = 1.0;
113: Gprev = 0;
114: Xprev = 0;
116: scale = 0;
118: H0 = 0;
119: H0default = TAO_TRUE;
120: }
124: TaoLMVMMat::~TaoLMVMMat()
125: {
126: int info;
128: info = S[0]->DestroyVecs(lm+1, S);
129: info = Y[0]->DestroyVecs(lm+1, Y);
130: info = TaoVecDestroy(D);
131: info = TaoVecDestroy(U);
132: info = TaoVecDestroy(V);
133: info = TaoVecDestroy(W);
134: info = TaoVecDestroy(P);
135: info = TaoVecDestroy(Q);
137: delete[] rho;
138: delete[] beta;
139: if (info) rho=0;
141: delete[] yy_history;
142: delete[] ys_history;
143: delete[] ss_history;
145: delete[] yy_rhistory;
146: delete[] ys_rhistory;
147: delete[] ss_rhistory;
148: }
152: int TaoLMVMMat::Refine(TaoLMVMMat *coarse, TaoMat *op, TaoVec *fineX, TaoVec *fineG)
153: {
154: double rhotemp, rhotol;
155: double y0temp, s0temp;
156: int info;
157: TaoInt i;
159: TaoFunctionBegin;
161: info = Reset(); CHKERRQ(info);
163: for (i = 0; i < coarse->lmnow; ++i) {
164: // Refine S[i] and Y[i]
165: info = op->Multiply(coarse->S[i], S[lmnow]); CHKERRQ(info);
166: info = op->Multiply(coarse->Y[i], Y[lmnow]); CHKERRQ(info);
168: // Check to see if refined vectors are fine
169: info = Y[lmnow]->Dot(S[lmnow], &rhotemp); CHKERRQ(info);
170: info = Y[lmnow]->Dot(Y[lmnow], &y0temp); CHKERRQ(info);
172: rhotol = eps * y0temp;
173: if (rhotemp > rhotol) {
174: rho[lmnow] = 1.0 / rhotemp;
176: // Add information to the scaling
177: switch(scaleType) {
178: case Scale_None:
179: break;
180:
181: case Scale_Scalar:
182: // Compute s^T s
183: info = S[lmnow]->Dot(S[lmnow], &s0temp); CHKERRQ(info);
185: // Save information for scalar scaling
186: yy_history[updates % scalar_history] = y0temp;
187: ys_history[updates % scalar_history] = rhotemp;
188: ss_history[updates % scalar_history] = s0temp;
190: sigma = coarse->sigma;
191: break;
193: case Scale_Broyden:
194: switch(rScaleType) {
195: case Rescale_None:
196: break;
197:
198: case Rescale_Scalar:
199: case Rescale_GL:
200: // Compute s^T s
201: info = S[lmnow]->Dot(S[lmnow], &s0temp); CHKERRQ(info);
203: // Save information for special cases of scalar rescaling
204: yy_rhistory[updates % rescale_history] = y0temp;
205: ys_rhistory[updates % rescale_history] = rhotemp;
206: ss_rhistory[updates % rescale_history] = s0temp;
207: break;
208: }
210: // Obtain diagonal scaling and ensure positive definite
211: info = op->Multiply(coarse->D, D); CHKERRQ(info);
212: info = D->AbsoluteValue(); CHKERRQ(info);
213: break;
214: }
216: // Successful update
217: ++updates;
218: ++lmnow;
219: }
220: else {
221: ++rejects;
222: }
223: }
225: // Save the number of iterations and previous values for x and g
226: // Need to use the actual value for the gradient here and not
227: // the refined version of the coarse gradient in order for the
228: // Wolfe conditions to guarantee the update is acceptable.
229: iter = coarse->iter;
230: info = Xprev->CopyFrom(fineX); CHKERRQ(info);
231: info = Gprev->CopyFrom(fineG); CHKERRQ(info);
232: TaoFunctionReturn(0);
233: }
237: int TaoLMVMMat::Reset()
238: {
239: TaoInt i;
241: TaoFunctionBegin;
242: Gprev = Y[lm];
243: Xprev = S[lm];
244: for (i = 0; i < lm; ++i) {
245: rho[i] = 0.0;
246: }
247: rho[0] = 1.0;
249: // Set the scaling and diagonal scaling matrix
250: switch(scaleType) {
251: case Scale_None:
252: sigma = 1.0;
253: break;
255: case Scale_Scalar:
256: sigma = delta;
257: break;
259: case Scale_Broyden:
260: D->SetToConstant(delta);
261: break;
262: }
264: iter = 0;
265: updates = 0;
266: lmnow = 0;
267: TaoFunctionReturn(0);
268: }
272: int TaoLMVMMat::Presolve()
273: {
274: TaoFunctionBegin;
275: TaoFunctionReturn(0);
276: }
280: int TaoLMVMMat::Solve(TaoVec *tb, TaoVec *dx, TaoTruth *tt)
281: {
282: double sq, yq, dd;
283: int info;
284: TaoInt ll;
285: TaoTruth scaled;
287: TaoFunctionBegin;
288: if (lmnow < 1) {
289: rho[0] = 1.0;
290: }
292: info = dx->CopyFrom(tb); CHKERRQ(info);
293: for (ll = 0; ll < lmnow; ++ll) {
294: info = dx->Dot(S[ll], &sq); CHKERRQ(info);
295: beta[ll] = sq * rho[ll];
296: info = dx->Axpy(-beta[ll], Y[ll]); CHKERRQ(info);
297: }
299: scaled = TAO_FALSE;
300: if (!scaled && H0) {
301: info = H0->Solve(dx, U, tt); CHKERRQ(info);
302: if (*tt) {
303: info = dx->Dot(U, &dd); CHKERRQ(info);
304: if ((dd > 0.0) && !TaoInfOrNaN(dd)) {
305: // Accept Hessian solve
306: info = dx->CopyFrom(U); CHKERRQ(info);
307: scaled = TAO_TRUE;
308: }
309: }
310: }
312: if (!scaled && scale) {
313: info = U->PointwiseMultiply(dx, scale); CHKERRQ(info);
314: info = dx->Dot(U, &dd); CHKERRQ(info);
315: if ((dd > 0.0) && !TaoInfOrNaN(dd)) {
316: // Accept scaling
317: info = dx->CopyFrom(U); CHKERRQ(info);
318: scaled = TAO_TRUE;
319: }
320: }
321:
322: if (!scaled) {
323: switch(scaleType) {
324: case Scale_None:
325: break;
327: case Scale_Scalar:
328: info = dx->Scale(sigma); CHKERRQ(info);
329: break;
330:
331: case Scale_Broyden:
332: info = dx->PointwiseMultiply(dx, D); CHKERRQ(info);
333: break;
334: }
335: }
337: for (ll = lmnow-1; ll >= 0; --ll) {
338: info = dx->Dot(Y[ll], &yq); CHKERRQ(info);
339: info = dx->Axpy(beta[ll] - yq * rho[ll], S[ll]); CHKERRQ(info);
340: }
342: *tt = TAO_TRUE;
343: TaoFunctionReturn(0);
344: }
348: int TaoLMVMMat::Update(TaoVec *x, TaoVec *g)
349: {
350: double rhotemp, rhotol;
351: double y0temp, s0temp;
352: double yDy, yDs, sDs;
353: double sigmanew, denom;
354: int info;
355: TaoInt i;
357: double yy_sum, ys_sum, ss_sum;
359: TaoFunctionBegin;
361: if (0 == iter) {
362: info = Reset(); CHKERRQ(info);
363: }
364: else {
365: info = Gprev->Aypx(-1.0, g); CHKERRQ(info);
366: info = Xprev->Aypx(-1.0, x); CHKERRQ(info);
368: info = Gprev->Dot(Xprev, &rhotemp); CHKERRQ(info);
369: info = Gprev->Dot(Gprev, &y0temp); CHKERRQ(info);
371: rhotol = eps * y0temp;
372: if (rhotemp > rhotol) {
373: ++updates;
375: lmnow = TaoMin(lmnow+1, lm);
376: for (i = lm-1; i >= 0; --i) {
377: S[i+1] = S[i];
378: Y[i+1] = Y[i];
379: rho[i+1] = rho[i];
380: }
381: S[0] = Xprev;
382: Y[0] = Gprev;
383: rho[0] = 1.0 / rhotemp;
385: // Compute the scaling
386: switch(scaleType) {
387: case Scale_None:
388: break;
390: case Scale_Scalar:
391: // Compute s^T s
392: info = Xprev->Dot(Xprev, &s0temp); CHKERRQ(info);
394: // Scalar is positive; safeguards are not required.
396: // Save information for scalar scaling
397: yy_history[(updates - 1) % scalar_history] = y0temp;
398: ys_history[(updates - 1) % scalar_history] = rhotemp;
399: ss_history[(updates - 1) % scalar_history] = s0temp;
401: // Compute summations for scalar scaling
402: yy_sum = 0; // No safeguard required; y^T y > 0
403: ys_sum = 0; // No safeguard required; y^T s > 0
404: ss_sum = 0; // No safeguard required; s^T s > 0
405: for (i = 0; i < TaoMin(updates, scalar_history); ++i) {
406: yy_sum += yy_history[i];
407: ys_sum += ys_history[i];
408: ss_sum += ss_history[i];
409: }
411: if (0.0 == s_alpha) {
412: // Safeguard ys_sum
413: if (0.0 == ys_sum) {
414: ys_sum = TAO_ZER_SAFEGUARD;
415: }
417: sigmanew = ss_sum / ys_sum;
418: }
419: else if (0.5 == s_alpha) {
420: // Safeguard yy_sum
421: if (0.0 == yy_sum) {
422: yy_sum = TAO_ZER_SAFEGUARD;
423: }
425: sigmanew = sqrt(ss_sum / yy_sum);
426: }
427: else if (1.0 == s_alpha) {
428: // Safeguard yy_sum
429: if (0.0 == yy_sum) {
430: yy_sum = TAO_ZER_SAFEGUARD;
431: }
433: sigmanew = ys_sum / yy_sum;
434: }
435: else {
436: denom = 2*s_alpha*yy_sum;
438: // Safeguard denom
439: if (0.0 == denom) {
440: denom = TAO_ZER_SAFEGUARD;
441: }
443: sigmanew = ((2*s_alpha-1)*ys_sum +
444: sqrt((2*s_alpha-1)*(2*s_alpha-1)*ys_sum*ys_sum -
445: 4*s_alpha*(s_alpha-1)*yy_sum*ss_sum)) / denom;
446: }
448: switch(limitType) {
449: case Limit_Average:
450: if (1.0 == mu) {
451: sigma = sigmanew;
452: }
453: else if (mu) {
454: sigma = mu * sigmanew + (1.0 - mu) * sigma;
455: }
456: break;
458: case Limit_Relative:
459: if (mu) {
460: sigma = TaoMid((1.0 - mu) * sigma, sigmanew, (1.0 + mu) * sigma);
461: }
462: break;
464: case Limit_Absolute:
465: if (nu) {
466: sigma = TaoMid(sigma - nu, sigmanew, sigma + nu);
467: }
468: break;
470: default:
471: sigma = sigmanew;
472: break;
473: }
474: break;
476: case Scale_Broyden:
477: // Original version
478: // Combine DFP and BFGS
480: // This code appears to be numerically unstable. We use the
481: // original version because this was used to generate all of
482: // the data and because it may be the least unstable of the
483: // bunch.
485: // P = Q = inv(D);
486: info = P->CopyFrom(D); CHKERRQ(info);
487: info = P->Reciprocal(); CHKERRQ(info);
488: info = Q->CopyFrom(P); CHKERRQ(info);
490: // V = y*y
491: info = V->PointwiseMultiply(Gprev, Gprev); CHKERRQ(info);
493: // W = inv(D)*s
494: info = W->PointwiseMultiply(Xprev, P); CHKERRQ(info);
495: info = W->Dot(Xprev, &sDs); CHKERRQ(info);
497: // Safeguard rhotemp and sDs
498: if (0.0 == rhotemp) {
499: rhotemp = TAO_ZER_SAFEGUARD;
500: }
502: if (0.0 == sDs) {
503: sDs = TAO_ZER_SAFEGUARD;
504: }
506: if (1.0 != phi) {
507: // BFGS portion of the update
508: // U = (inv(D)*s)*(inv(D)*s)
509: info = U->PointwiseMultiply(W, W); CHKERRQ(info);
511: // Assemble
512: info = P->Axpy(1.0 / rhotemp, V); CHKERRQ(info);
513: info = P->Axpy(-1.0 / sDs, U); CHKERRQ(info);
514: }
516: if (0.0 != phi) {
517: // DFP portion of the update
518: // U = inv(D)*s*y
519: info = U->PointwiseMultiply(W, Gprev); CHKERRQ(info);
521: // Assemble
522: info = Q->Axpy(1.0 / rhotemp + sDs / (rhotemp * rhotemp), V); CHKERRQ(info);
523: info = Q->Axpy(-2.0 / rhotemp, U); CHKERRQ(info);
524: }
526: if (0.0 == phi) {
527: info = U->CopyFrom(P); CHKERRQ(info);
528: }
529: else if (1.0 == phi) {
530: info = U->CopyFrom(Q); CHKERRQ(info);
531: }
532: else {
533: // Broyden update
534: info = U->Waxpby(1.0 - phi, P, phi, Q); CHKERRQ(info);
535: }
537: // Obtain inverse and ensure positive definite
538: info = U->Reciprocal(); CHKERRQ(info);
539: info = U->AbsoluteValue(); CHKERRQ(info);
541: // Checking the diagonal scaling for not a number and infinity
542: // should not be necessary for the Broyden update
543: // info = U->Dot(U, &sDs); CHKERRQ(info);
544: // if (sDs != sDs) {
545: // // not a number
546: // info = U->SetToConstant(TAO_ZER_SAFEGUARD); CHKERRQ(info);
547: // }
548: // else if ((sDs - sDs) != 0.0)) {
549: // // infinity
550: // info = U->SetToConstant(TAO_INF_SAFEGUARD); CHKERRQ(info);
551: // }
553: switch(rScaleType) {
554: case Rescale_None:
555: break;
557: case Rescale_Scalar:
558: case Rescale_GL:
559: if (rScaleType == Rescale_GL) {
560: // Gilbert and Lemarachal use the old diagonal
561: info = P->CopyFrom(D); CHKERRQ(info);
562: }
563: else {
564: // The default version uses the current diagonal
565: info = P->CopyFrom(U); CHKERRQ(info);
566: }
568: // Compute s^T s
569: info = Xprev->Dot(Xprev, &s0temp); CHKERRQ(info);
571: // Save information for special cases of scalar rescaling
572: yy_rhistory[(updates - 1) % rescale_history] = y0temp;
573: ys_rhistory[(updates - 1) % rescale_history] = rhotemp;
574: ss_rhistory[(updates - 1) % rescale_history] = s0temp;
576: if (0.0 == r_beta) {
577: if (1 == TaoMin(updates, rescale_history)) {
578: // Compute summations for scalar scaling
579: info = W->PointwiseDivide(S[0], P); CHKERRQ(info);
580:
581: info = W->Dot(Y[0], &ys_sum); CHKERRQ(info);
582: info = W->Dot(W, &ss_sum); CHKERRQ(info);
583:
584: yy_sum += yy_rhistory[0];
585: }
586: else {
587: info = Q->CopyFrom(P); CHKERRQ(info);
588: info = Q->Reciprocal(); CHKERRQ(info);
590: // Compute summations for scalar scaling
591: yy_sum = 0; // No safeguard required
592: ys_sum = 0; // No safeguard required
593: ss_sum = 0; // No safeguard required
594: for (i = 0; i < TaoMin(updates, rescale_history); ++i) {
595: info = W->PointwiseMultiply(S[i], Q); CHKERRQ(info);
596: info = W->Dot(Y[i], &yDs); CHKERRQ(info);
597: ys_sum += yDs;
599: info = W->Dot(W, &sDs); CHKERRQ(info);
600: ss_sum += sDs;
601:
602: yy_sum += yy_rhistory[i];
603: }
604: }
605: }
606: else if (0.5 == r_beta) {
607: if (1 == TaoMin(updates, rescale_history)) {
608: info = V->PointwiseMultiply(Y[0], P); CHKERRQ(info);
609: info = V->Dot(Y[0], &yy_sum); CHKERRQ(info);
611: info = W->PointwiseDivide(S[0], P); CHKERRQ(info);
612: info = W->Dot(S[0], &ss_sum); CHKERRQ(info);
614: ys_sum = ys_rhistory[0];
615: }
616: else {
617: info = Q->CopyFrom(P); CHKERRQ(info);
618: info = Q->Reciprocal(); CHKERRQ(info);
620: // Compute summations for scalar scaling
621: yy_sum = 0; // No safeguard required
622: ys_sum = 0; // No safeguard required
623: ss_sum = 0; // No safeguard required
624: for (i = 0; i < TaoMin(updates, rescale_history); ++i) {
625: info = V->PointwiseMultiply(Y[i], P); CHKERRQ(info);
626: info = V->Dot(Y[i], &yDy); CHKERRQ(info);
627: yy_sum += yDy;
629: info = W->PointwiseMultiply(S[i], Q); CHKERRQ(info);
630: info = W->Dot(S[i], &sDs); CHKERRQ(info);
631: ss_sum += sDs;
633: ys_sum += ys_rhistory[i];
634: }
635: }
636: }
637: else if (1.0 == r_beta) {
638: // Compute summations for scalar scaling
639: yy_sum = 0; // No safeguard required
640: ys_sum = 0; // No safeguard required
641: ss_sum = 0; // No safeguard required
642: for (i = 0; i < TaoMin(updates, rescale_history); ++i) {
643: info = V->PointwiseMultiply(Y[i], P); CHKERRQ(info);
644: info = V->Dot(S[i], &yDs); CHKERRQ(info);
645: ys_sum += yDs;
647: info = V->Dot(V, &yDy); CHKERRQ(info);
648: yy_sum += yDy;
650: ss_sum += ss_rhistory[i];
651: }
652: }
653: else {
654: info = Q->CopyFrom(P); CHKERRQ(info);
656: info = P->Pow(r_beta); CHKERRQ(info);
657: info = Q->PointwiseDivide(P, Q); CHKERRQ(info);
659: // Compute summations for scalar scaling
660: yy_sum = 0; // No safeguard required
661: ys_sum = 0; // No safeguard required
662: ss_sum = 0; // No safeguard required
663: for (i = 0; i < TaoMin(updates, rescale_history); ++i) {
664: info = V->PointwiseMultiply(P, Y[i]); CHKERRQ(info);
665: info = W->PointwiseMultiply(Q, S[i]); CHKERRQ(info);
667: info = V->Dot(V, &yDy); CHKERRQ(info);
668: info = V->Dot(W, &yDs); CHKERRQ(info);
669: info = W->Dot(W, &sDs); CHKERRQ(info);
671: yy_sum += yDy;
672: ys_sum += yDs;
673: ss_sum += sDs;
674: }
675: }
677: if (0.0 == r_alpha) {
678: // Safeguard ys_sum
679: if (0.0 == ys_sum) {
680: ys_sum = TAO_ZER_SAFEGUARD;
681: }
683: sigmanew = ss_sum / ys_sum;
684: }
685: else if (0.5 == r_alpha) {
686: // Safeguard yy_sum
687: if (0.0 == yy_sum) {
688: yy_sum = TAO_ZER_SAFEGUARD;
689: }
691: sigmanew = sqrt(ss_sum / yy_sum);
692: }
693: else if (1.0 == r_alpha) {
694: // Safeguard yy_sum
695: if (0.0 == yy_sum) {
696: yy_sum = TAO_ZER_SAFEGUARD;
697: }
699: sigmanew = ys_sum / yy_sum;
700: }
701: else {
702: denom = 2*r_alpha*yy_sum;
704: // Safeguard denom
705: if (0.0 == denom) {
706: denom = TAO_ZER_SAFEGUARD;
707: }
709: sigmanew = ((2*r_alpha-1)*ys_sum +
710: sqrt((2*r_alpha-1)*(2*r_alpha-1)*ys_sum*ys_sum -
711: 4*r_alpha*(r_alpha-1)*yy_sum*ss_sum)) / denom;
712: }
714: // If Q has small values, then Q^(r_beta - 1)
715: // can have very large values. Hence, ys_sum
716: // and ss_sum can be infinity. In this case,
717: // sigmanew can either be not-a-number or infinity.
719: if (TaoInfOrNaN(sigmanew)) {
720: // sigmanew is not-a-number; skip rescaling
721: }
722: else if (!sigmanew) {
723: // sigmanew is zero; this is a bad case; skip rescaling
724: }
725: else {
726: // sigmanew is positive
727: info = U->Scale(sigmanew); CHKERRQ(info);
728: }
729: break;
730: }
732: // Modify for previous information
733: switch(limitType) {
734: case Limit_Average:
735: if (1.0 == mu) {
736: info = D->CopyFrom(U); CHKERRQ(info);
737: }
738: else if (mu) {
739: info = D->Axpby(mu, U, 1.0 - mu); CHKERRQ(info);
740: }
741: break;
742:
743: case Limit_Relative:
744: if (mu) {
745: info = P->ScaleCopyFrom(1.0 - mu, D); CHKERRQ(info);
746: info = Q->ScaleCopyFrom(1.0 + mu, D); CHKERRQ(info);
747: info = D->Median(P, U, Q); CHKERRQ(info);
748: }
749: break;
751: case Limit_Absolute:
752: if (nu) {
753: info = P->CopyFrom(D); CHKERRQ(info);
754: info = P->AddConstant(-nu); CHKERRQ(info);
755: info = Q->CopyFrom(D); CHKERRQ(info);
756: info = Q->AddConstant(nu); CHKERRQ(info);
757: info = D->Median(P, U, Q); CHKERRQ(info);
758: }
759: break;
761: default:
762: info = D->CopyFrom(U); CHKERRQ(info);
763: break;
764: }
765: break;
766: }
768: Xprev = S[lm]; Gprev = Y[lm];
769: }
770: else {
771: ++rejects;
772: }
773: }
774:
775: ++iter;
776: info = Xprev->CopyFrom(x); CHKERRQ(info);
777: info = Gprev->CopyFrom(g); CHKERRQ(info);
778: TaoFunctionReturn(0);
779: }
783: int TaoLMVMMat::View()
784: {
785: TaoFunctionBegin;
786: TaoFunctionReturn(0);
787: }
791: int TaoLMVMMat::SetDelta(double d)
792: {
793: TaoFunctionBegin;
794: delta = TaoAbsDouble(d);
795: delta = TaoMax(delta_min, delta);
796: delta = TaoMin(delta_max, delta);
797: TaoFunctionReturn(0);
798: }
802: int TaoLMVMMat::SetScale(TaoVec *s)
803: {
804: TaoFunctionBegin;
805: scale = s;
806: TaoFunctionReturn(0);
807: }
811: int TaoLMVMMat::SetH0(TaoMat *HH0)
812: {
813: TaoFunctionBegin;
814: if (H0) {
815: H0 = 0;
816: }
817: H0default = TAO_TRUE;
819: if (HH0) {
820: H0 = HH0;
821: H0default = TAO_FALSE;
822: }
823: TaoFunctionReturn(0);
824: }
828: int TaoLMVMMat::GetX0(TaoVec *x)
829: {
830: int i,info;
832: TaoFunctionBegin;
833: info = x->CopyFrom(Xprev); CHKERRQ(info);
834: for (i = 0; i < lmnow; ++i) {
835: info = x->Axpy(-1.0, S[i]); CHKERRQ(info);
836: }
837: TaoFunctionReturn(0);
838: }
842: int TaoLMVMMat::InitialApproximation(TaoVec *x)
843: {
844: TaoFunctionBegin;
845: TaoFunctionReturn(0);
846: }