[occt.git] / src / FEmTool / FEmTool_ProfileMatrix.cxx

// Created on: 1997-10-31
// Created by: Roman BORISOV
// Copyright (c) 1997-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#define No_Standard_RangeError
#define No_Standard_OutOfRange
#define No_Standard_DimensionError

#include <FEmTool_ProfileMatrix.ixx>
#include <gp.hxx>


//=======================================================================
//function : :FEmTool_ProfileMatrix
//purpose  :
//=======================================================================
FEmTool_ProfileMatrix::FEmTool_ProfileMatrix(const TColStd_Array1OfInteger& FirstIndexes) 
                       : profile(1, 2, 1, FirstIndexes.Length())
{
  Standard_Integer i, j, k, l;
  profile(1, 1) = 0;
  profile(2, 1) = 1;
  for(i = 2; i <= FirstIndexes.Length(); i++) {
    profile(1, i) = i - FirstIndexes(i);
    profile(2, i) = profile(2, i-1) + profile(1, i) + 1;
  }
  NextCoeff = new TColStd_HArray1OfInteger(1, profile(2, FirstIndexes.Length()));
 
  for(i = 1, k = 1; i <= FirstIndexes.Length(); i++)
    for(j = FirstIndexes(i); j <= i; j++) {
      for(l = i+1; l <= FirstIndexes.Length() && j < FirstIndexes(l); l++);

      if(l > FirstIndexes.Length()) NextCoeff->SetValue(k, 0);
      else NextCoeff->SetValue(k, l);
      k++;
    }
  

  ProfileMatrix = new TColStd_HArray1OfReal(1, profile(2, FirstIndexes.Length()));
  SMatrix = new TColStd_HArray1OfReal(1, profile(2, FirstIndexes.Length()));
  IsDecomp = Standard_False;
}

//=======================================================================
//function : Init
//purpose  :
//=======================================================================

 void FEmTool_ProfileMatrix::Init(const Standard_Real Value) 
{
  ProfileMatrix->Init(Value);
  IsDecomp = Standard_False;
}

//=======================================================================
//function : ChangeValue
//purpose  : 
//=======================================================================

 Standard_Real& FEmTool_ProfileMatrix::ChangeValue(const Standard_Integer I,
						   const Standard_Integer J)
{
  Standard_Integer Ind;
  Ind = I-J;
  if (Ind < 0) {
    Ind = -Ind;
    Standard_OutOfRange_Raise_if(Ind>profile(1, J),
				  "FEmTool_ProfileMatrix::ChangeValue");
    Ind = profile(2, J) - Ind;
  }
  else {
    Standard_OutOfRange_Raise_if(Ind>profile(1, I),
				  "FEmTool_ProfileMatrix::ChangeValue");
    Ind = profile(2, I)-Ind;
  }
  return  ProfileMatrix->ChangeValue(Ind);
}

//=======================================================================
//function : Decompose
//purpose  : Choleski's decomposition
//=======================================================================
 Standard_Boolean FEmTool_ProfileMatrix::Decompose() 
{
  Standard_Integer i, j, k, ik, jk, DiagAddr, CurrAddr, Kmin, Kj;
  Standard_Real Sum, a, Eps = 1.e-32;

  SMatrix->Init(0.);
  Standard_Real * SMA =  &SMatrix->ChangeValue(1);
  SMA--;
  const Standard_Real * PM = &ProfileMatrix->Value(1);
  PM--;

  for(j = 1; j <= RowNumber(); j++) {
    DiagAddr = profile(2, j);
    Kj = j - profile(1, j);
    Sum = 0;
    for(k = DiagAddr - profile(1, j); k < DiagAddr; k++)
      Sum += SMA[k] * SMA[k];

    a = PM[DiagAddr] - Sum;
    if(a < Eps) {
      return IsDecomp = Standard_False;// Matrix is not positive defined
    }
    a = Sqrt(a);
    SMA[DiagAddr] = a;
  
    CurrAddr = DiagAddr;
    while ((i = NextCoeff->Value(CurrAddr)) > 0) {
      CurrAddr = profile(2, i) - (i - j);
    
      // Computation of Sum of S  .S  for k = 1,..,j-1
      //                        ik  jk
      Sum = 0;
      Kmin = Max((i - profile(1, i)), Kj);
      ik = profile(2, i) - i + Kmin;
      jk = DiagAddr - j + Kmin;
      for(k = Kmin; k <j; k++,ik++,jk++) {
	Sum += SMA[ik]*SMA[jk];
      }  
      SMA[CurrAddr] = (PM[CurrAddr] - Sum)/a;
    }
  }
  return IsDecomp = Standard_True;
}

//=======================================================================
//function : Solve
//purpose  : Resolution of the system S*t(S)X = B
//=======================================================================
 void FEmTool_ProfileMatrix::Solve(const math_Vector& B,math_Vector& X) const
{
  if (!IsDecomp) StdFail_NotDone::Raise("Decomposition must be done");
 
  Standard_Integer i, j, jj,DiagAddr, CurrAddr;
  Standard_Real Sum;

  Standard_Real * x = &X(X.Lower());
  x--;
  Standard_Real * b = &B(B.Lower());
  b--;
  const Standard_Real * SMA = &SMatrix->Value(1);
  SMA --;
  const Standard_Integer * NC = &NextCoeff->Value(1);
  NC--;

// Resolution of Sw = B;
  for(i = 1; i <= RowNumber(); i++) {
    DiagAddr = profile(2, i);
    Sum = 0;
    for(j = i - profile(1, i), jj = DiagAddr - (i - j); 
	j < i; j++, jj++) 
      Sum += SMA[jj]* x[j];
    x[i] = (b[i] - Sum)/SMA[DiagAddr];
  }

// Resolution of t(S)X = w;
  for(i = ColNumber(); i >= 1; i--) {
    DiagAddr = profile(2, i);
    j = NC[DiagAddr];
    Sum = 0;
    while(j > 0) {
      CurrAddr = profile(2, j) - (j-i);
      Sum += SMA[CurrAddr]*x[j];
      j = NC[CurrAddr];
    }
   x[i] = (x[i] - Sum)/SMA[DiagAddr]; 
  }
}

 Standard_Boolean FEmTool_ProfileMatrix::Prepare() 
{
  Standard_NotImplemented::Raise("FEmTool_ProfileMatrix::Prepare");
  return Standard_False;  
}

// void FEmTool_ProfileMatrix::Solve(const math_Vector& B,const math_Vector& Init,math_Vector& X,math_Vector& Residual,const Standard_Real Tolerance,const Standard_Integer NbIterations) const
 void FEmTool_ProfileMatrix::Solve(const math_Vector& ,const math_Vector& ,math_Vector& ,math_Vector& ,const Standard_Real ,const Standard_Integer ) const
{
  Standard_NotImplemented::Raise("FEmTool_ProfileMatrix::Solve");
}


//=======================================================================
//function : Multiplied
//purpose  : MX = H*X
//=======================================================================
 void FEmTool_ProfileMatrix::Multiplied(const math_Vector& X,math_Vector& MX) const
{
  Standard_Integer i, j, jj, DiagAddr, CurrAddr;
  Standard_Real * m = &MX(MX.Lower());
  m--;
  Standard_Real * x = &X(X.Lower());
  x--;
  const Standard_Real * PM = &ProfileMatrix->Value(1);
  PM--;
  const Standard_Integer * NC = &NextCoeff->Value(1);
  NC--;

  for(i = 1; i <= RowNumber(); i++) {
    DiagAddr = profile(2, i);
    m[i] = 0;
    for(j = i - profile(1, i), jj = DiagAddr - (i - j); 
	j <= i; j++, jj++)
      m[i] += PM[jj]*x[j];

    CurrAddr = DiagAddr;
    for(j = NC[CurrAddr]; j > 0; j = NC[CurrAddr]) {
      CurrAddr = profile(2, j) - (j-i);
      m[i] += PM[CurrAddr]*x[j];
    }
  }
}

 Standard_Integer FEmTool_ProfileMatrix::RowNumber() const
{
  return profile.RowLength();
} 

 Standard_Integer FEmTool_ProfileMatrix::ColNumber() const
{
  return profile.RowLength();
}

Standard_Boolean FEmTool_ProfileMatrix::IsInProfile(const Standard_Integer i, 
						    const Standard_Integer j) const
{
  if (j <= i) {
    if ((i - j) <= profile(1, i)) return Standard_True;
    else return Standard_False;
  }
  else if ((j - i) <= profile(1, j)) return Standard_True;

  return Standard_False;
}

 void FEmTool_ProfileMatrix::OutM() const
{
  Standard_Integer i, j;
  cout<<"Matrix A"<<endl;
  for(i = 1; i <= RowNumber(); i++) {
    for(j = 1; j < i - profile(1, i); j++) 
      cout<<"0 ";

    for(j = profile(2, i) - profile(1, i); j <= profile(2, i); j++) 
      cout<<ProfileMatrix->Value(j)<<" ";
    cout<<endl;
  }

  cout<<"NextCoeff"<<endl;
  for(i = 1; i <= profile(2, RowNumber()); i++)
    cout<<NextCoeff->Value(i)<<" ";
  cout<<endl;
}

 void FEmTool_ProfileMatrix::OutS() const
{
  Standard_Integer i, j;
  cout<<"Matrix S"<<endl;
  for(i = 1; i <= RowNumber(); i++) {
    for(j = 1; j < i - profile(1, i); j++) 
      cout<<"0 ";

    for(j = profile(2, i) - profile(1, i); j <= profile(2, i); j++) 
      cout<<SMatrix->Value(j)<<" ";
    cout<<endl;
  }
}
Commit	Line	Data
b311480e	1	// Created on: 1997-10-31
	2	// Created by: Roman BORISOV
	3	// Copyright (c) 1997-1999 Matra Datavision
973c2be1	4	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	5	//
973c2be1	6	// This file is part of Open CASCADE Technology software library.
b311480e	7	//
d5f74e42	8	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	9	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	10	// by the Free Software Foundation, with special exception defined in the file
	11	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	12	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	13	//
973c2be1	14	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	15	// commercial license or contractual agreement.
7fd59977	16
	17	#define No_Standard_RangeError
	18	#define No_Standard_OutOfRange
	19	#define No_Standard_DimensionError
	20
	21	#include <FEmTool_ProfileMatrix.ixx>
	22	#include <gp.hxx>
	23
	24
	25	//=======================================================================
	26	//function : :FEmTool_ProfileMatrix
	27	//purpose :
	28	//=======================================================================
b311480e	29	FEmTool_ProfileMatrix::FEmTool_ProfileMatrix(const TColStd_Array1OfInteger& FirstIndexes)
7fd59977	30	: profile(1, 2, 1, FirstIndexes.Length())
	31	{
	32	Standard_Integer i, j, k, l;
	33	profile(1, 1) = 0;
	34	profile(2, 1) = 1;
	35	for(i = 2; i <= FirstIndexes.Length(); i++) {
	36	profile(1, i) = i - FirstIndexes(i);
	37	profile(2, i) = profile(2, i-1) + profile(1, i) + 1;
	38	}
	39	NextCoeff = new TColStd_HArray1OfInteger(1, profile(2, FirstIndexes.Length()));
	40
	41	for(i = 1, k = 1; i <= FirstIndexes.Length(); i++)
	42	for(j = FirstIndexes(i); j <= i; j++) {
	43	for(l = i+1; l <= FirstIndexes.Length() && j < FirstIndexes(l); l++);
	44
	45	if(l > FirstIndexes.Length()) NextCoeff->SetValue(k, 0);
	46	else NextCoeff->SetValue(k, l);
	47	k++;
	48	}
	49
	50
	51	ProfileMatrix = new TColStd_HArray1OfReal(1, profile(2, FirstIndexes.Length()));
	52	SMatrix = new TColStd_HArray1OfReal(1, profile(2, FirstIndexes.Length()));
	53	IsDecomp = Standard_False;
	54	}
	55
	56	//=======================================================================
	57	//function : Init
	58	//purpose :
	59	//=======================================================================
	60
	61	void FEmTool_ProfileMatrix::Init(const Standard_Real Value)
	62	{
	63	ProfileMatrix->Init(Value);
	64	IsDecomp = Standard_False;
	65	}
	66
	67	//=======================================================================
	68	//function : ChangeValue
	69	//purpose :
	70	//=======================================================================
	71
	72	Standard_Real& FEmTool_ProfileMatrix::ChangeValue(const Standard_Integer I,
	73	const Standard_Integer J)
	74	{
	75	Standard_Integer Ind;
	76	Ind = I-J;
	77	if (Ind < 0) {
	78	Ind = -Ind;
	79	Standard_OutOfRange_Raise_if(Ind>profile(1, J),
	80	"FEmTool_ProfileMatrix::ChangeValue");
	81	Ind = profile(2, J) - Ind;
	82	}
	83	else {
	84	Standard_OutOfRange_Raise_if(Ind>profile(1, I),
	85	"FEmTool_ProfileMatrix::ChangeValue");
	86	Ind = profile(2, I)-Ind;
	87	}
	88	return ProfileMatrix->ChangeValue(Ind);
	89	}
	90
	91	//=======================================================================
	92	//function : Decompose
	93	//purpose : Choleski's decomposition
94	//=======================================================================
95	Standard_Boolean FEmTool_ProfileMatrix::Decompose()
96	{
97	Standard_Integer i, j, k, ik, jk, DiagAddr, CurrAddr, Kmin, Kj;
98	Standard_Real Sum, a, Eps = 1.e-32;
99
100	SMatrix->Init(0.);
101	Standard_Real * SMA = &SMatrix->ChangeValue(1);
102	SMA--;
103	const Standard_Real * PM = &ProfileMatrix->Value(1);
104	PM--;
105
106	for(j = 1; j <= RowNumber(); j++) {
107	DiagAddr = profile(2, j);
108	Kj = j - profile(1, j);
109	Sum = 0;
110	for(k = DiagAddr - profile(1, j); k < DiagAddr; k++)
111	Sum += SMA[k] * SMA[k];
112
113	a = PM[DiagAddr] - Sum;
114	if(a < Eps) {
115	return IsDecomp = Standard_False;// Matrix is not positive defined
116	}
117	a = Sqrt(a);
118	SMA[DiagAddr] = a;
119
120	CurrAddr = DiagAddr;
121	while ((i = NextCoeff->Value(CurrAddr)) > 0) {
122	CurrAddr = profile(2, i) - (i - j);
123
124	// Computation of Sum of S .S for k = 1,..,j-1
125	// ik jk
126	Sum = 0;
127	Kmin = Max((i - profile(1, i)), Kj);
128	ik = profile(2, i) - i + Kmin;
129	jk = DiagAddr - j + Kmin;
130	for(k = Kmin; k <j; k++,ik++,jk++) {
131	Sum += SMA[ik]*SMA[jk];
132	}
133	SMA[CurrAddr] = (PM[CurrAddr] - Sum)/a;
134	}
135	}
136	return IsDecomp = Standard_True;
137	}
138
139	//=======================================================================
140	//function : Solve
141	//purpose : Resolution of the system S*t(S)X = B
142	//=======================================================================
143	void FEmTool_ProfileMatrix::Solve(const math_Vector& B,math_Vector& X) const
144	{
145	if (!IsDecomp) StdFail_NotDone::Raise("Decomposition must be done");
146
147	Standard_Integer i, j, jj,DiagAddr, CurrAddr;
148	Standard_Real Sum;
149
150	Standard_Real * x = &X(X.Lower());
151	x--;
152	Standard_Real * b = &B(B.Lower());
153	b--;
154	const Standard_Real * SMA = &SMatrix->Value(1);
155	SMA --;
156	const Standard_Integer * NC = &NextCoeff->Value(1);
157	NC--;
158
159	// Resolution of Sw = B;
160	for(i = 1; i <= RowNumber(); i++) {
161	DiagAddr = profile(2, i);
162	Sum = 0;
163	for(j = i - profile(1, i), jj = DiagAddr - (i - j);
164	j < i; j++, jj++)
165	Sum += SMA[jj]* x[j];
166	x[i] = (b[i] - Sum)/SMA[DiagAddr];
167	}
168
169	// Resolution of t(S)X = w;
170	for(i = ColNumber(); i >= 1; i--) {
171	DiagAddr = profile(2, i);
172	j = NC[DiagAddr];
173	Sum = 0;
174	while(j > 0) {
175	CurrAddr = profile(2, j) - (j-i);
176	Sum += SMA[CurrAddr]*x[j];
177	j = NC[CurrAddr];
178	}
179	x[i] = (x[i] - Sum)/SMA[DiagAddr];
180	}
181	}
182
183	Standard_Boolean FEmTool_ProfileMatrix::Prepare()
184	{
185	Standard_NotImplemented::Raise("FEmTool_ProfileMatrix::Prepare");
186	return Standard_False;
187	}
188
189	// void FEmTool_ProfileMatrix::Solve(const math_Vector& B,const math_Vector& Init,math_Vector& X,math_Vector& Residual,const Standard_Real Tolerance,const Standard_Integer NbIterations) const
190	void FEmTool_ProfileMatrix::Solve(const math_Vector& ,const math_Vector& ,math_Vector& ,math_Vector& ,const Standard_Real ,const Standard_Integer ) const
191	{
192	Standard_NotImplemented::Raise("FEmTool_ProfileMatrix::Solve");
193	}
194
195
196	//=======================================================================
197	//function : Multiplied
198	//purpose : MX = H*X
199	//=======================================================================
200	void FEmTool_ProfileMatrix::Multiplied(const math_Vector& X,math_Vector& MX) const
201	{
202	Standard_Integer i, j, jj, DiagAddr, CurrAddr;
203	Standard_Real * m = &MX(MX.Lower());
204	m--;
205	Standard_Real * x = &X(X.Lower());
206	x--;
207	const Standard_Real * PM = &ProfileMatrix->Value(1);
208	PM--;
209	const Standard_Integer * NC = &NextCoeff->Value(1);
210	NC--;
211
212	for(i = 1; i <= RowNumber(); i++) {
213	DiagAddr = profile(2, i);
214	m[i] = 0;
215	for(j = i - profile(1, i), jj = DiagAddr - (i - j);
216	j <= i; j++, jj++)
217	m[i] += PM[jj]*x[j];
218
219	CurrAddr = DiagAddr;
220	for(j = NC[CurrAddr]; j > 0; j = NC[CurrAddr]) {
221	CurrAddr = profile(2, j) - (j-i);
222	m[i] += PM[CurrAddr]*x[j];
223	}
224	}
225	}
226
227	Standard_Integer FEmTool_ProfileMatrix::RowNumber() const
228	{
229	return profile.RowLength();
230	}
231
232	Standard_Integer FEmTool_ProfileMatrix::ColNumber() const
233	{
234	return profile.RowLength();
235	}
236
237	Standard_Boolean FEmTool_ProfileMatrix::IsInProfile(const Standard_Integer i,
238	const Standard_Integer j) const
239	{
240	if (j <= i) {
241	if ((i - j) <= profile(1, i)) return Standard_True;
242	else return Standard_False;
243	}
244	else if ((j - i) <= profile(1, j)) return Standard_True;
245
246	return Standard_False;
247	}
248
249	void FEmTool_ProfileMatrix::OutM() const
250	{
251	Standard_Integer i, j;
252	cout<<"Matrix A"<<endl;
253	for(i = 1; i <= RowNumber(); i++) {
254	for(j = 1; j < i - profile(1, i); j++)
255	cout<<"0 ";
256
257	for(j = profile(2, i) - profile(1, i); j <= profile(2, i); j++)
258	cout<<ProfileMatrix->Value(j)<<" ";
259	cout<<endl;
260	}
261
262	cout<<"NextCoeff"<<endl;
263	for(i = 1; i <= profile(2, RowNumber()); i++)
264	cout<<NextCoeff->Value(i)<<" ";
265	cout<<endl;
266	}
267
268	void FEmTool_ProfileMatrix::OutS() const
269	{
270	Standard_Integer i, j;
271	cout<<"Matrix S"<<endl;
272	for(i = 1; i <= RowNumber(); i++) {
273	for(j = 1; j < i - profile(1, i); j++)
274	cout<<"0 ";
275
276	for(j = profile(2, i) - profile(1, i); j <= profile(2, i); j++)
277	cout<<SMatrix->Value(j)<<" ";
278	cout<<endl;
279	}
280	}