Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/trlinsolve.cs @ 2806

Visit:

Last change on this file since 2806 was 2806, checked in by gkronber, 15 years ago
Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 39.0 KB

Rev	Line
[2806]	1	/*************************************************************************
	2	This file is a part of ALGLIB project.
	3
	4	>>> SOURCE LICENSE >>>
	5	This program is free software; you can redistribute it and/or modify
	6	it under the terms of the GNU General Public License as published by
	7	the Free Software Foundation (www.fsf.org); either version 2 of the
	8	License, or (at your option) any later version.
	9
	10	This program is distributed in the hope that it will be useful,
	11	but WITHOUT ANY WARRANTY; without even the implied warranty of
	12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	13	GNU General Public License for more details.
	14
	15	A copy of the GNU General Public License is available at
	16	http://www.fsf.org/licensing/licenses
	17
	18	>>> END OF LICENSE >>>
	19	*************************************************************************/
	20
	21	using System;
	22
	23	namespace alglib
	24	{
	25	public class trlinsolve
	26	{
	27	/*************************************************************************
	28	Utility subroutine performing the "safe" solution of system of linear
	29	equations with triangular coefficient matrices.
	30
	31	The subroutine uses scaling and solves the scaled system Ax=sb (where s
	32	is a scalar value) instead of A*x=b, choosing s so that x can be
	33	represented by a floating-point number. The closer the system gets to a
	34	singular, the less s is. If the system is singular, s=0 and x contains the
	35	non-trivial solution of equation A*x=0.
	36
	37	The feature of an algorithm is that it could not cause an overflow or a
	38	division by zero regardless of the matrix used as the input.
	39
	40	The algorithm can solve systems of equations with upper/lower triangular
	41	matrices, with/without unit diagonal, and systems of type Ax=b or A'x=b
	42	(where A' is a transposed matrix A).
	43
	44	Input parameters:
	45	A - system matrix. Array whose indexes range within [0..N-1, 0..N-1].
	46	N - size of matrix A.
	47	X - right-hand member of a system.
	48	Array whose index ranges within [0..N-1].
	49	IsUpper - matrix type. If it is True, the system matrix is the upper
	50	triangular and is located in the corresponding part of
	51	matrix A.
	52	Trans - problem type. If it is True, the problem to be solved is
	53	A'x=b, otherwise it is Ax=b.
	54	Isunit - matrix type. If it is True, the system matrix has a unit
	55	diagonal (the elements on the main diagonal are not used
	56	in the calculation process), otherwise the matrix is considered
	57	to be a general triangular matrix.
	58
	59	Output parameters:
	60	X - solution. Array whose index ranges within [0..N-1].
	61	S - scaling factor.
	62
	63	-- LAPACK auxiliary routine (version 3.0) --
	64	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	65	Courant Institute, Argonne National Lab, and Rice University
	66	June 30, 1992
	67	*************************************************************************/
	68	public static void rmatrixtrsafesolve(ref double[,] a,
	69	int n,
	70	ref double[] x,
	71	ref double s,
	72	bool isupper,
	73	bool istrans,
	74	bool isunit)
	75	{
	76	bool normin = new bool();
	77	double[] cnorm = new double[0];
	78	double[,] a1 = new double[0,0];
	79	double[] x1 = new double[0];
	80	int i = 0;
	81	int i_ = 0;
	82	int i1_ = 0;
	83
	84
	85	//
	86	// From 0-based to 1-based
	87	//
	88	normin = false;
	89	a1 = new double[n+1, n+1];
	90	x1 = new double[n+1];
	91	for(i=1; i<=n; i++)
	92	{
	93	i1_ = (0) - (1);
	94	for(i_=1; i_<=n;i_++)
	95	{
	96	a1[i,i_] = a[i-1,i_+i1_];
	97	}
	98	}
	99	i1_ = (0) - (1);
	100	for(i_=1; i_<=n;i_++)
	101	{
	102	x1[i_] = x[i_+i1_];
	103	}
	104
	105	//
	106	// Solve 1-based
	107	//
	108	safesolvetriangular(ref a1, n, ref x1, ref s, isupper, istrans, isunit, normin, ref cnorm);
	109
	110	//
	111	// From 1-based to 0-based
	112	//
	113	i1_ = (1) - (0);
	114	for(i_=0; i_<=n-1;i_++)
	115	{
	116	x[i_] = x1[i_+i1_];
	117	}
	118	}
	119
	120
	121	/*************************************************************************
	122	Obsolete 1-based subroutine.
	123	See RMatrixTRSafeSolve for 0-based replacement.
	124	*************************************************************************/
	125	public static void safesolvetriangular(ref double[,] a,
	126	int n,
	127	ref double[] x,
	128	ref double s,
	129	bool isupper,
	130	bool istrans,
	131	bool isunit,
	132	bool normin,
	133	ref double[] cnorm)
	134	{
	135	int i = 0;
	136	int imax = 0;
	137	int j = 0;
	138	int jfirst = 0;
	139	int jinc = 0;
	140	int jlast = 0;
	141	int jm1 = 0;
	142	int jp1 = 0;
	143	int ip1 = 0;
	144	int im1 = 0;
	145	int k = 0;
	146	int flg = 0;
	147	double v = 0;
	148	double vd = 0;
	149	double bignum = 0;
	150	double grow = 0;
	151	double rec = 0;
	152	double smlnum = 0;
	153	double sumj = 0;
	154	double tjj = 0;
	155	double tjjs = 0;
	156	double tmax = 0;
	157	double tscal = 0;
	158	double uscal = 0;
	159	double xbnd = 0;
	160	double xj = 0;
	161	double xmax = 0;
	162	bool notran = new bool();
	163	bool upper = new bool();
	164	bool nounit = new bool();
	165	int i_ = 0;
	166
	167	upper = isupper;
	168	notran = !istrans;
	169	nounit = !isunit;
	170
	171	//
	172	// Quick return if possible
	173	//
	174	if( n==0 )
	175	{
	176	return;
	177	}
	178
	179	//
	180	// Determine machine dependent parameters to control overflow.
	181	//
	182	smlnum = AP.Math.MinRealNumber/(AP.Math.MachineEpsilon*2);
	183	bignum = 1/smlnum;
	184	s = 1;
	185	if( !normin )
	186	{
	187	cnorm = new double[n+1];
	188
	189	//
	190	// Compute the 1-norm of each column, not including the diagonal.
	191	//
	192	if( upper )
	193	{
	194
	195	//
	196	// A is upper triangular.
	197	//
	198	for(j=1; j<=n; j++)
	199	{
	200	v = 0;
	201	for(k=1; k<=j-1; k++)
	202	{
	203	v = v+Math.Abs(a[k,j]);
	204	}
	205	cnorm[j] = v;
	206	}
	207	}
	208	else
	209	{
	210
	211	//
	212	// A is lower triangular.
	213	//
	214	for(j=1; j<=n-1; j++)
	215	{
	216	v = 0;
	217	for(k=j+1; k<=n; k++)
	218	{
	219	v = v+Math.Abs(a[k,j]);
	220	}
	221	cnorm[j] = v;
	222	}
	223	cnorm[n] = 0;
	224	}
	225	}
	226
	227	//
	228	// Scale the column norms by TSCAL if the maximum element in CNORM is
	229	// greater than BIGNUM.
	230	//
	231	imax = 1;
	232	for(k=2; k<=n; k++)
	233	{
	234	if( (double)(cnorm[k])>(double)(cnorm[imax]) )
	235	{
	236	imax = k;
	237	}
	238	}
	239	tmax = cnorm[imax];
	240	if( (double)(tmax)<=(double)(bignum) )
	241	{
	242	tscal = 1;
	243	}
	244	else
	245	{
	246	tscal = 1/(smlnum*tmax);
	247	for(i_=1; i_<=n;i_++)
	248	{
	249	cnorm[i_] = tscal*cnorm[i_];
	250	}
	251	}
	252
	253	//
	254	// Compute a bound on the computed solution vector to see if the
	255	// Level 2 BLAS routine DTRSV can be used.
	256	//
	257	j = 1;
	258	for(k=2; k<=n; k++)
	259	{
	260	if( (double)(Math.Abs(x[k]))>(double)(Math.Abs(x[j])) )
	261	{
	262	j = k;
	263	}
	264	}
	265	xmax = Math.Abs(x[j]);
	266	xbnd = xmax;
	267	if( notran )
	268	{
	269
	270	//
	271	// Compute the growth in A * x = b.
	272	//
	273	if( upper )
	274	{
	275	jfirst = n;
	276	jlast = 1;
	277	jinc = -1;
	278	}
	279	else
	280	{
	281	jfirst = 1;
	282	jlast = n;
	283	jinc = 1;
	284	}
	285	if( (double)(tscal)!=(double)(1) )
	286	{
	287	grow = 0;
	288	}
	289	else
	290	{
	291	if( nounit )
	292	{
	293
	294	//
	295	// A is non-unit triangular.
	296	//
	297	// Compute GROW = 1/G(j) and XBND = 1/M(j).
	298	// Initially, G(0) = max{x(i), i=1,...,n}.
	299	//
	300	grow = 1/Math.Max(xbnd, smlnum);
	301	xbnd = grow;
	302	j = jfirst;
	303	while( jinc>0 & j<=jlast \| jinc<0 & j>=jlast )
	304	{
	305
	306	//
	307	// Exit the loop if the growth factor is too small.
	308	//
	309	if( (double)(grow)<=(double)(smlnum) )
	310	{
	311	break;
	312	}
	313
	314	//
	315	// M(j) = G(j-1) / abs(A(j,j))
	316	//
	317	tjj = Math.Abs(a[j,j]);
	318	xbnd = Math.Min(xbnd, Math.Min(1, tjj)*grow);
	319	if( (double)(tjj+cnorm[j])>=(double)(smlnum) )
	320	{
	321
	322	//
	323	// G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
	324	//
	325	grow = grow*(tjj/(tjj+cnorm[j]));
	326	}
	327	else
	328	{
	329
	330	//
	331	// G(j) could overflow, set GROW to 0.
	332	//
	333	grow = 0;
	334	}
	335	if( j==jlast )
	336	{
	337	grow = xbnd;
	338	}
	339	j = j+jinc;
	340	}
	341	}
	342	else
	343	{
	344
	345	//
	346	// A is unit triangular.
	347	//
	348	// Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
	349	//
	350	grow = Math.Min(1, 1/Math.Max(xbnd, smlnum));
	351	j = jfirst;
	352	while( jinc>0 & j<=jlast \| jinc<0 & j>=jlast )
	353	{
	354
	355	//
	356	// Exit the loop if the growth factor is too small.
	357	//
	358	if( (double)(grow)<=(double)(smlnum) )
	359	{
	360	break;
	361	}
	362
	363	//
	364	// G(j) = G(j-1)*( 1 + CNORM(j) )
	365	//
	366	grow = grow*(1/(1+cnorm[j]));
	367	j = j+jinc;
	368	}
	369	}
	370	}
	371	}
	372	else
	373	{
	374
	375	//
	376	// Compute the growth in A' * x = b.
	377	//
	378	if( upper )
	379	{
	380	jfirst = 1;
	381	jlast = n;
	382	jinc = 1;
	383	}
	384	else
	385	{
	386	jfirst = n;
	387	jlast = 1;
	388	jinc = -1;
	389	}
	390	if( (double)(tscal)!=(double)(1) )
	391	{
	392	grow = 0;
	393	}
	394	else
	395	{
	396	if( nounit )
	397	{
	398
	399	//
	400	// A is non-unit triangular.
	401	//
	402	// Compute GROW = 1/G(j) and XBND = 1/M(j).
	403	// Initially, M(0) = max{x(i), i=1,...,n}.
	404	//
	405	grow = 1/Math.Max(xbnd, smlnum);
	406	xbnd = grow;
	407	j = jfirst;
	408	while( jinc>0 & j<=jlast \| jinc<0 & j>=jlast )
	409	{
	410
	411	//
	412	// Exit the loop if the growth factor is too small.
	413	//
	414	if( (double)(grow)<=(double)(smlnum) )
	415	{
	416	break;
	417	}
	418
	419	//
	420	// G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
	421	//
	422	xj = 1+cnorm[j];
	423	grow = Math.Min(grow, xbnd/xj);
	424
	425	//
	426	// M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
	427	//
	428	tjj = Math.Abs(a[j,j]);
	429	if( (double)(xj)>(double)(tjj) )
	430	{
	431	xbnd = xbnd*(tjj/xj);
	432	}
	433	if( j==jlast )
	434	{
	435	grow = Math.Min(grow, xbnd);
	436	}
	437	j = j+jinc;
	438	}
	439	}
	440	else
	441	{
	442
	443	//
	444	// A is unit triangular.
	445	//
	446	// Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
	447	//
	448	grow = Math.Min(1, 1/Math.Max(xbnd, smlnum));
	449	j = jfirst;
	450	while( jinc>0 & j<=jlast \| jinc<0 & j>=jlast )
	451	{
	452
	453	//
	454	// Exit the loop if the growth factor is too small.
	455	//
	456	if( (double)(grow)<=(double)(smlnum) )
	457	{
	458	break;
	459	}
	460
	461	//
	462	// G(j) = ( 1 + CNORM(j) )*G(j-1)
	463	//
	464	xj = 1+cnorm[j];
	465	grow = grow/xj;
	466	j = j+jinc;
	467	}
	468	}
	469	}
	470	}
	471	if( (double)(grow*tscal)>(double)(smlnum) )
	472	{
	473
	474	//
	475	// Use the Level 2 BLAS solve if the reciprocal of the bound on
	476	// elements of X is not too small.
	477	//
	478	if( upper & notran \| !upper & !notran )
	479	{
	480	if( nounit )
	481	{
	482	vd = a[n,n];
	483	}
	484	else
	485	{
	486	vd = 1;
	487	}
	488	x[n] = x[n]/vd;
	489	for(i=n-1; i>=1; i--)
	490	{
	491	ip1 = i+1;
	492	if( upper )
	493	{
	494	v = 0.0;
	495	for(i_=ip1; i_<=n;i_++)
	496	{
	497	v += a[i,i_]*x[i_];
	498	}
	499	}
	500	else
	501	{
	502	v = 0.0;
	503	for(i_=ip1; i_<=n;i_++)
	504	{
	505	v += a[i_,i]*x[i_];
	506	}
	507	}
	508	if( nounit )
	509	{
	510	vd = a[i,i];
	511	}
	512	else
	513	{
	514	vd = 1;
	515	}
	516	x[i] = (x[i]-v)/vd;
	517	}
	518	}
	519	else
	520	{
	521	if( nounit )
	522	{
	523	vd = a[1,1];
	524	}
	525	else
	526	{
	527	vd = 1;
	528	}
	529	x[1] = x[1]/vd;
	530	for(i=2; i<=n; i++)
	531	{
	532	im1 = i-1;
	533	if( upper )
	534	{
	535	v = 0.0;
	536	for(i_=1; i_<=im1;i_++)
	537	{
	538	v += a[i_,i]*x[i_];
	539	}
	540	}
	541	else
	542	{
	543	v = 0.0;
	544	for(i_=1; i_<=im1;i_++)
	545	{
	546	v += a[i,i_]*x[i_];
	547	}
	548	}
	549	if( nounit )
	550	{
	551	vd = a[i,i];
	552	}
	553	else
	554	{
	555	vd = 1;
	556	}
	557	x[i] = (x[i]-v)/vd;
	558	}
	559	}
	560	}
	561	else
	562	{
	563
	564	//
	565	// Use a Level 1 BLAS solve, scaling intermediate results.
	566	//
	567	if( (double)(xmax)>(double)(bignum) )
	568	{
	569
	570	//
	571	// Scale X so that its components are less than or equal to
	572	// BIGNUM in absolute value.
	573	//
	574	s = bignum/xmax;
	575	for(i_=1; i_<=n;i_++)
	576	{
	577	x[i_] = s*x[i_];
	578	}
	579	xmax = bignum;
	580	}
	581	if( notran )
	582	{
	583
	584	//
	585	// Solve A * x = b
	586	//
	587	j = jfirst;
	588	while( jinc>0 & j<=jlast \| jinc<0 & j>=jlast )
	589	{
	590
	591	//
	592	// Compute x(j) = b(j) / A(j,j), scaling x if necessary.
	593	//
	594	xj = Math.Abs(x[j]);
	595	flg = 0;
	596	if( nounit )
	597	{
	598	tjjs = a[j,j]*tscal;
	599	}
	600	else
	601	{
	602	tjjs = tscal;
	603	if( (double)(tscal)==(double)(1) )
	604	{
	605	flg = 100;
	606	}
	607	}
	608	if( flg!=100 )
	609	{
	610	tjj = Math.Abs(tjjs);
	611	if( (double)(tjj)>(double)(smlnum) )
	612	{
	613
	614	//
	615	// abs(A(j,j)) > SMLNUM:
	616	//
	617	if( (double)(tjj)<(double)(1) )
	618	{
	619	if( (double)(xj)>(double)(tjj*bignum) )
	620	{
	621
	622	//
	623	// Scale x by 1/b(j).
	624	//
	625	rec = 1/xj;
	626	for(i_=1; i_<=n;i_++)
	627	{
	628	x[i_] = rec*x[i_];
	629	}
	630	s = s*rec;
	631	xmax = xmax*rec;
	632	}
	633	}
	634	x[j] = x[j]/tjjs;
	635	xj = Math.Abs(x[j]);
	636	}
	637	else
	638	{
	639	if( (double)(tjj)>(double)(0) )
	640	{
	641
	642	//
	643	// 0 < abs(A(j,j)) <= SMLNUM:
	644	//
	645	if( (double)(xj)>(double)(tjj*bignum) )
	646	{
	647
	648	//
	649	// Scale x by (1/abs(x(j)))abs(A(j,j))BIGNUM
	650	// to avoid overflow when dividing by A(j,j).
	651	//
	652	rec = tjj*bignum/xj;
	653	if( (double)(cnorm[j])>(double)(1) )
	654	{
	655
	656	//
	657	// Scale by 1/CNORM(j) to avoid overflow when
	658	// multiplying x(j) times column j.
	659	//
	660	rec = rec/cnorm[j];
	661	}
	662	for(i_=1; i_<=n;i_++)
	663	{
	664	x[i_] = rec*x[i_];
	665	}
	666	s = s*rec;
	667	xmax = xmax*rec;
	668	}
	669	x[j] = x[j]/tjjs;
	670	xj = Math.Abs(x[j]);
	671	}
	672	else
	673	{
	674
	675	//
	676	// A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
	677	// scale = 0, and compute a solution to A*x = 0.
	678	//
	679	for(i=1; i<=n; i++)
	680	{
	681	x[i] = 0;
	682	}
	683	x[j] = 1;
	684	xj = 1;
	685	s = 0;
	686	xmax = 0;
	687	}
	688	}
	689	}
	690
	691	//
	692	// Scale x if necessary to avoid overflow when adding a
	693	// multiple of column j of A.
	694	//
	695	if( (double)(xj)>(double)(1) )
	696	{
	697	rec = 1/xj;
	698	if( (double)(cnorm[j])>(double)((bignum-xmax)*rec) )
	699	{
	700
	701	//
	702	// Scale x by 1/(2*abs(x(j))).
	703	//
	704	rec = rec*0.5;
	705	for(i_=1; i_<=n;i_++)
	706	{
	707	x[i_] = rec*x[i_];
	708	}
	709	s = s*rec;
	710	}
	711	}
	712	else
	713	{
	714	if( (double)(xj*cnorm[j])>(double)(bignum-xmax) )
	715	{
	716
	717	//
	718	// Scale x by 1/2.
	719	//
	720	for(i_=1; i_<=n;i_++)
	721	{
	722	x[i_] = 0.5*x[i_];
	723	}
	724	s = s*0.5;
	725	}
	726	}
	727	if( upper )
	728	{
	729	if( j>1 )
	730	{
	731
	732	//
	733	// Compute the update
	734	// x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
	735	//
	736	v = x[j]*tscal;
	737	jm1 = j-1;
	738	for(i_=1; i_<=jm1;i_++)
	739	{
	740	x[i_] = x[i_] - v*a[i_,j];
	741	}
	742	i = 1;
	743	for(k=2; k<=j-1; k++)
	744	{
	745	if( (double)(Math.Abs(x[k]))>(double)(Math.Abs(x[i])) )
	746	{
	747	i = k;
	748	}
	749	}
	750	xmax = Math.Abs(x[i]);
	751	}
	752	}
	753	else
	754	{
	755	if( j<n )
	756	{
	757
	758	//
	759	// Compute the update
	760	// x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
	761	//
	762	jp1 = j+1;
	763	v = x[j]*tscal;
	764	for(i_=jp1; i_<=n;i_++)
	765	{
	766	x[i_] = x[i_] - v*a[i_,j];
	767	}
	768	i = j+1;
	769	for(k=j+2; k<=n; k++)
	770	{
	771	if( (double)(Math.Abs(x[k]))>(double)(Math.Abs(x[i])) )
	772	{
	773	i = k;
	774	}
	775	}
	776	xmax = Math.Abs(x[i]);
	777	}
	778	}
	779	j = j+jinc;
	780	}
	781	}
	782	else
	783	{
	784
	785	//
	786	// Solve A' * x = b
	787	//
	788	j = jfirst;
	789	while( jinc>0 & j<=jlast \| jinc<0 & j>=jlast )
	790	{
	791
	792	//
	793	// Compute x(j) = b(j) - sum A(k,j)*x(k).
	794	// k<>j
	795	//
	796	xj = Math.Abs(x[j]);
	797	uscal = tscal;
	798	rec = 1/Math.Max(xmax, 1);
	799	if( (double)(cnorm[j])>(double)((bignum-xj)*rec) )
	800	{
	801
	802	//
	803	// If x(j) could overflow, scale x by 1/(2*XMAX).
	804	//
	805	rec = rec*0.5;
	806	if( nounit )
	807	{
	808	tjjs = a[j,j]*tscal;
	809	}
	810	else
	811	{
	812	tjjs = tscal;
	813	}
	814	tjj = Math.Abs(tjjs);
	815	if( (double)(tjj)>(double)(1) )
	816	{
	817
	818	//
	819	// Divide by A(j,j) when scaling x if A(j,j) > 1.
	820	//
	821	rec = Math.Min(1, rec*tjj);
	822	uscal = uscal/tjjs;
	823	}
	824	if( (double)(rec)<(double)(1) )
	825	{
	826	for(i_=1; i_<=n;i_++)
	827	{
	828	x[i_] = rec*x[i_];
	829	}
	830	s = s*rec;
	831	xmax = xmax*rec;
	832	}
	833	}
	834	sumj = 0;
	835	if( (double)(uscal)==(double)(1) )
	836	{
	837
	838	//
	839	// If the scaling needed for A in the dot product is 1,
	840	// call DDOT to perform the dot product.
	841	//
	842	if( upper )
	843	{
	844	if( j>1 )
	845	{
	846	jm1 = j-1;
	847	sumj = 0.0;
	848	for(i_=1; i_<=jm1;i_++)
	849	{
	850	sumj += a[i_,j]*x[i_];
	851	}
	852	}
	853	else
	854	{
	855	sumj = 0;
	856	}
	857	}
	858	else
	859	{
	860	if( j<n )
	861	{
	862	jp1 = j+1;
	863	sumj = 0.0;
	864	for(i_=jp1; i_<=n;i_++)
	865	{
	866	sumj += a[i_,j]*x[i_];
	867	}
	868	}
	869	}
	870	}
	871	else
	872	{
	873
	874	//
	875	// Otherwise, use in-line code for the dot product.
	876	//
	877	if( upper )
	878	{
	879	for(i=1; i<=j-1; i++)
	880	{
	881	v = a[i,j]*uscal;
	882	sumj = sumj+v*x[i];
	883	}
	884	}
	885	else
	886	{
	887	if( j<n )
	888	{
	889	for(i=j+1; i<=n; i++)
	890	{
	891	v = a[i,j]*uscal;
	892	sumj = sumj+v*x[i];
	893	}
	894	}
	895	}
	896	}
	897	if( (double)(uscal)==(double)(tscal) )
	898	{
	899
	900	//
	901	// Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
	902	// was not used to scale the dotproduct.
	903	//
	904	x[j] = x[j]-sumj;
	905	xj = Math.Abs(x[j]);
	906	flg = 0;
	907	if( nounit )
	908	{
	909	tjjs = a[j,j]*tscal;
	910	}
	911	else
	912	{
	913	tjjs = tscal;
	914	if( (double)(tscal)==(double)(1) )
	915	{
	916	flg = 150;
	917	}
	918	}
	919
	920	//
	921	// Compute x(j) = x(j) / A(j,j), scaling if necessary.
	922	//
	923	if( flg!=150 )
	924	{
	925	tjj = Math.Abs(tjjs);
	926	if( (double)(tjj)>(double)(smlnum) )
	927	{
	928
	929	//
	930	// abs(A(j,j)) > SMLNUM:
	931	//
	932	if( (double)(tjj)<(double)(1) )
	933	{
	934	if( (double)(xj)>(double)(tjj*bignum) )
	935	{
	936
	937	//
	938	// Scale X by 1/abs(x(j)).
	939	//
	940	rec = 1/xj;
	941	for(i_=1; i_<=n;i_++)
	942	{
	943	x[i_] = rec*x[i_];
	944	}
	945	s = s*rec;
	946	xmax = xmax*rec;
	947	}
	948	}
	949	x[j] = x[j]/tjjs;
	950	}
	951	else
	952	{
	953	if( (double)(tjj)>(double)(0) )
	954	{
	955
	956	//
	957	// 0 < abs(A(j,j)) <= SMLNUM:
	958	//
	959	if( (double)(xj)>(double)(tjj*bignum) )
	960	{
	961
	962	//
	963	// Scale x by (1/abs(x(j)))abs(A(j,j))BIGNUM.
	964	//
	965	rec = tjj*bignum/xj;
	966	for(i_=1; i_<=n;i_++)
	967	{
	968	x[i_] = rec*x[i_];
	969	}
	970	s = s*rec;
	971	xmax = xmax*rec;
	972	}
	973	x[j] = x[j]/tjjs;
	974	}
	975	else
	976	{
	977
	978	//
	979	// A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
	980	// scale = 0, and compute a solution to A'*x = 0.
	981	//
	982	for(i=1; i<=n; i++)
	983	{
	984	x[i] = 0;
	985	}
	986	x[j] = 1;
	987	s = 0;
	988	xmax = 0;
	989	}
	990	}
	991	}
	992	}
	993	else
	994	{
	995
	996	//
	997	// Compute x(j) := x(j) / A(j,j) - sumj if the dot
	998	// product has already been divided by 1/A(j,j).
	999	//
	1000	x[j] = x[j]/tjjs-sumj;
	1001	}
	1002	xmax = Math.Max(xmax, Math.Abs(x[j]));
	1003	j = j+jinc;
	1004	}
	1005	}
	1006	s = s/tscal;
	1007	}
	1008
	1009	//
	1010	// Scale the column norms by 1/TSCAL for return.
	1011	//
	1012	if( (double)(tscal)!=(double)(1) )
	1013	{
	1014	v = 1/tscal;
	1015	for(i_=1; i_<=n;i_++)
	1016	{
	1017	cnorm[i_] = v*cnorm[i_];
	1018	}
	1019	}
	1020	}
	1021	}
	1022	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences