001 package numal.highprecision;
002
003 import java.util.*;
004 import numal.highprecision.*;
005
006 public class Analytic_problems extends Object {
007
008 public static void marquardt(int m, int n, double par[],
009 double g[], double v[][],
010 AP_marquardt_methods method,
011 double in[], double out[])
012 {
013 boolean emergency;
014 int maxfe,fe,it,i,j,err;
015 double vv,ww,w,mu,res,fpar,fparpres,lambda,lambdamin,p,pw,
016 reltolres,abstolres,temp;
017 double em[] = new double[8];
018 double val[] = new double[n+1];
019 double b[] = new double[n+1];
020 double bb[] = new double[n+1];
021 double parpres[] = new double[n+1];
022 double jac[][] = new double[m+1][n+1];
023
024 lambda=0.0;
025 vv=10.0;
026 w=0.5;
027 mu=0.01;
028 ww = (in[6] < 1.0e-7) ? 1.0e-8 : 1.0e-1*in[6];
029 em[0]=em[2]=em[6]=in[0];
030 em[4]=10*n;
031 reltolres=in[3];
032 abstolres=in[4]*in[4];
033 maxfe=(int)in[5];
034 err=0;
035 fe=it=1;
036 p=fpar=res=0.0;
037 pw = -Math.log(ww*in[0])/2.30;
038 if (!method.funct(m,n,par,g)) {
039 err=3;
040 out[4]=fe;
041 out[5]=it-1;
042 out[1]=err;
043 return;
044 }
045 fpar=Basic.vecvec(1,m,0,g,g);
046 out[3]=Math.sqrt(fpar);
047 emergency=false;
048 it=1;
049 do {
050 method.jacobian(m,n,par,g,jac);
051 i=Linear_algebra.qrisngvaldec(jac,m,n,val,v,em);
052 if (it == 1)
053 lambda=in[6]*Basic.vecvec(1,n,0,val,val);
054 else
055 if (p == 0.0) lambda *= w;
056 for (i=1; i<=n; i++) b[i]=val[i]*Basic.tamvec(1,m,i,jac,g);
057 while (true) {
058 for (i=1; i<=n; i++) bb[i]=b[i]/(val[i]*val[i]+lambda);
059 for (i=1; i<=n; i++)
060 parpres[i]=par[i]-Basic.matvec(1,n,i,v,bb);
061 fe++;
062 if (fe >= maxfe)
063 err=1;
064 else
065 if (!method.funct(m,n,parpres,g)) err=2;
066 if (err != 0) {
067 emergency=true;
068 break;
069 }
070 fparpres=Basic.vecvec(1,m,0,g,g);
071 res=fpar-fparpres;
072 if (res < mu*Basic.vecvec(1,n,0,b,bb)) {
073 p += 1.0;
074 lambda *= vv;
075 if (p == 1.0) {
076 lambdamin=ww*Basic.vecvec(1,n,0,val,val);
077 if (lambda < lambdamin) lambda=lambdamin;
078 }
079 if (p >= pw) {
080 err=4;
081 emergency=true;
082 break;
083 }
084 } else {
085 Basic.dupvec(1,n,0,par,parpres);
086 fpar=fparpres;
087 break;
088 }
089 }
090 if (emergency) break;
091 it++;
092 } while (fpar > abstolres && res > reltolres*fpar+abstolres);
093 for (i=1; i<=n; i++)
094 Basic.mulcol(1,n,i,i,jac,v,1.0/(val[i]+in[0]));
095 for (i=1; i<=n; i++)
096 for (j=1; j<=i; j++)
097 v[i][j]=v[j][i]=Basic.mattam(1,n,i,j,jac,jac);
098 lambda=lambdamin=val[1];
099 for (i=2; i<=n; i++)
100 if (val[i] > lambda)
101 lambda=val[i];
102 else
103 if (val[i] < lambdamin) lambdamin=val[i];
104 temp=lambda/(lambdamin+in[0]);
105 out[7]=temp*temp;
106 out[2]=Math.sqrt(fpar);
107 out[6]=Math.sqrt(res+fpar)-out[2];
108 out[4]=fe;
109 out[5]=it-1;
110 out[1]=err;
111 }
112
113 }