Пример:
Код: Выделить всё
A =
9.0 0.0 0.0 3.0 1.0 0.0 1.0
0.0 11.0 2.0 1.0 0.0 0.0 2.0
0.0 2.0 10.0 2.0 0.0 0.0 0.0
3.0 1.0 2.0 9.0 1.0 0.0 0.0
1.0 0.0 0.0 1.0 12.0 0.0 1.0
0.0 0.0 0.0 0.0 0.0 8.0 0.0
1.0 2.0 0.0 0.0 1.0 0.0 8.0
------------- Разреженная матрица ------------
Формат хранения: CSIR
AN: 9.0 3.0 1.0 1.0 11.0 2.0 1.0 2.0 10.0 2.0 9.0 1.0 12.0 1.0 8.0 8.0
JA: 0.0 3.0 4.0 6.0 1.0 2.0 3.0 6.0 2.0 3.0 3.0 4.0 4.0 6.0 5.0 6.0
IA: 0.0 4.0 8.0 10.0 12.0 14.0 15.0 16.0
-----------------------------------------------
Код: Выделить всё
/**
* Неполная факторизация Холецкого разреженной матрицы
* @param AN - массив, содержащий все ненулевые над диагональные элементы матрицы A, включая и диагональные.
* Каждая новая строка этого массива начинается с диагонального элемента.
* @param JA - одномерный массив, который содержит столько же элементов, сколько и
* массив AN и для каждого из них указывает, в каком столбце находиться данный
* элемент.
* @param IA - одномерный массив, элементы которого указывают на позиции, с
* которых начинается описание очередной строки.
* @return u - массив, содержащий ненулевые наддиагональные елементы, получившиеся после неполной факторизации Холецкого
*/
public static double[] IC_Decomposition(double[] AN,int[] JA,int[] IA,double[] u ){
int i,j,k;
double sum2=0, m=0,r=0;
int nz =AN.length;//количество ненулевых элементов, не считая диагональные
int n = IA.length-1;//количество строк (столбцов), матрицы
double[] sum = new double[n];
u[0]=Math.sqrt(AN[0]);
//находим элементы 0 - ой строки
for (i=IA[0];i<IA[1];i++) {
u[i]=AN[i]/u[0];
sum[JA[i]]=u[i]*u[i];//запоминаем квадраты получившихся элементов
}
//Ищем остальные факторизованные элементы , начиная с 1-ой строки
for (i=1;i<n;i++){
u[IA[i]]=Math.sqrt(AN[IA[i]]-sum[i]);//находим диагональные элементы
for (j=IA[i]+1;j<IA[i+1];j++){
u[j] = AN[j];
//Поиск cуммы элементов a[k][i] и a[k][j].
/////////////////////////////////////////////////////////////////
for ( k = 0; k < i; k++)
{
int n1 = -1;
for (int ii = IA[k]+1; ii < IA[k + 1]; ii++)
{
if (JA[ii] == i) n1 = ii;
if (JA[ii] == JA[j])
{
if (n1 >= 0) u[j] -= u[n1]*u[ii];
break;
}
}
}
u[j]/=u[IA[i]];//Находим остальные над диагональные элементы
sum[JA[j]]+=u[j]*u[j];
}
}
return u;
}
Вот его код:
Код: Выделить всё
static void
ICFactor(const MV_Vector<int> &pntr, const MV_Vector<int> &indx,
MV_Vector<double> &val)
{
int d, g, h, i, j, k, n = pntr.size() - 1;
double z;
for (k = 0; k < n - 1; k++) {
d = pntr[k];
z = val[d] = sqrt(val[d]);
for (i = d + 1; i < pntr[k+1]; i++)
val[i] /= z;
for (i = d + 1; i < pntr[k+1]; i++) {
z = val[i];
h = indx[i];
g = i;
for (j = pntr[h] ; j < pntr[h+1]; j++)
for ( ; g < pntr[k+1] && indx[g+1] <= indx[j]; g++)
if (indx[g] == indx[j])
val[j] -= z * val[g];
}
}
d = pntr[n-1];
val[d] = sqrt(val[d]);
}