В данный момент занимаюсь написанием предобуславливателя к итерационоому решателю Conjuate Gradiuent (метод сопряженных градиентов) для sparse - матриц. Матрица храниться в 4-х одномерных массивах (разреженный строчный формат для симметричных матриц CSIR).
Пример:
Код:
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: 3.0 1.0 1.0 2.0 1.0 2.0 2.0 1.0 1.0
JA: 3.0 4.0 6.0 2.0 3.0 6.0 3.0 4.0 6.0
IA: 0.0 3.0 6.0 7.0 8.0 9.0 9.0 9.0
AD: 9.0 11.0 10.0 9.0 12.0 8.0 8.0
-----------------------------------------------
В качестве улутшателя взял Неполное разложение Холецкого (Incomplete Cholesky). Написал непосредственно само неполное разложение но я его реализовал так что оно , вследствии дополнительного поиска элементов выполняется очень долго. Хотелось бы ускорить процесс факторизации а то толку от предобуславливания ни какого не будет при большеразмерных матрицах.
Код:
/**
* Неполная факторизация Холецкого разреженной матрицы
* @param AN - массив, содержащий все ненулевые наддиагональные элементы матрицы A
* @param AD - массив, содержащий диагональные элементы матрицы A
* @param JA - одномерный массив, который содержит столько же элементов, сколько и
* массив AN и для каждого из них указывает, в каком столбце находиться данный элемент.
* @param IA - одномерный массив, элементы которого указывают на позиции, с которых *начинается описание очередной строки.
* @return u - массив, содержащий ненулевые наддиагональные елементы, получившиеся после неполной факторизации Холецкого
*/
public static double[] IC_Decomposition(double[] AN,double[] AD,int[] JA,int[] IA ){
int i,j,k;
double sum2=0, m=0,r=0;
int nz =AN.length;//количество ненулевых элементов, не считая диагональные
int n = IA.length-1;//количество строк (столбцов), матрицы
double[] u = new double[nz];
double[] udiag = new double[n];//массив, содержащий диагональные элементы матрицы A, получившиеся после неполной факторизации Холецкого
double[] sum = new double[n];
udiag[0]=Math.sqrt(AD[0]);
//находим элементы 0 - ой строки
for (i=IA[0];i<IA[1];i++) {
u[i]=AN[i]/udiag[0];
sum[JA[i]]=u[i]*u[i];//запоминаем квадраты получившихся элементов
}
//Ищем остальные факторизованные элементы , начиная с 1-ой строки
for (i=1;i<n;i++){
udiag[i]=Math.sqrt(AD[i]-sum[i]);//находим диагональные элементы
for (j=IA[i];j<IA[i+1];j++){
sum2=0;
//Поиск элементов a[k][i] и a[k][j]. Здесь на до бы придумать что-нибудь по проще. Слишком долго ищет!!!
/////////////////////////////////////////////////////////////////
for(k=0;k<i;k++) {
m=FindElement(k,i,u,IA,JA); // ищем элемент a[k][i]
r=FindElement(k,JA[j],u,IA,JA); // ищем элемент a[k][j]
if (m!=0 && r!=0) sum2+=m*r; //если оба не нулевые то перемножаем
}
///////////////////////////////////////////////////////////////
u[j]=(AN[j]-sum2)/udiag[i];//Находим остальные наддиагональные элементы
sum[JA[j]]+=u[j]*u[j];
}
}
return u;
}
/**
* Поиск элемента в разреженной матрице
* @param row - номер строки
* @param col - номер столбца
* @param u - массив в котором ищем
* @param IA -
* @param JA -
* @return elem - найденный элемент
*/
public static double FindElement(int row,int col,double[] u,int[] IA,int[] JA ){
int ii;
double elem=0.0;
for (ii=IA[row];ii<IA[row+1];ii++){
if (col==JA[ii]) {
elem = u[ii];
return elem;
}
}
return elem;
}
u =
1.0
0.3333333333333333
0.3333333333333333
0.6030226891555273
0.30151134457776363
0.6030226891555273
0.5857074126574228
0.24236755930688816
0.2584356424246762
udiag =
3.0
3.3166247903554
3.1042492870843406
2.750643149492325
3.439498052781032
2.8284271247461903
2.731018774006702
Ребята помогите придумать как найти ненулевые элементы a[k][i] и a[k][j] без открывания доп. циклов, или дайте пример на реализацию Неполного Холецкого для матриц , хранящихся в формате CSIR.
Нашел реализацию IC для sparse матриц в библиотеке ITL 4.0.0.1. Ту же самую реализацию нашел в библиотеке IML++.
Вот она:
Код:
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]);
}
Как я понял массив val это аналог моего AN и соответственно pntr - IA , indx - JA.
Только здесь матрица A храниться несколько в другом формате:
Пример:
Код:
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 в ITL
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
-----------------------------------------------
Я переписал этот алгоритм из библиотеки IML++ и у меня получилось вот что:
Код:
public static double[] IC_Decomposition2(double[] AN,int[] JA,int[] IA ){
int d, g, h, i, j, k;
double z;
int nz = AN.length;
int n = IA.length-1;
for (k = 0; k < n - 1; k++) {
d =(int) IA[k];
z = AN[d] = Math.sqrt(AN[d]);
for (i = d + 1; i < IA[k+1]; i++)
AN[i] /= z;
for (i = d + 1; i < IA[k+1]; i++) {
z = AN[i];
h = JA[i];
g = i;
for (j =IA[h] ; j < IA[h+1]; j++)
for ( ; g < IA[k+1] && JA[g+1] <= JA[j]; g++)
if (JA[g] == JA[j])
AN[j] -= z * AN[g];
}
}
d = IA[n-1];
AN[d] = Math.sqrt(AN[d]);
return AN;
}
AN (factorization)=
3.0
1.0
0.3333333333333333
0.3333333333333333
3.3166247903554
0.6030226891555273
0.30151134457776363
0.6030226891555273
3.1622776601683795
0.6324555320336759
2.932575659723036
0.34099716973523675
3.4472773213410837
0.2578524458667825
2.8284271247461903
2.7310738989293286
AN - массив ненулевых наддиагональных и диагональных элементов, получившихся после факторизации.
По скорости этот алгоритм намного быстрее моего(и по скорости он меня устраивает), но его решение не точное. Если вы посмотрите на решение , полученное в моём алгоритме, то увидите разницу.(напр. 0.5857(прав.) и 0.6324(неправ.) ; 2.75(прав.) и 2.93(неправ) и.т.д. )
Тут надо либо изменить поиск элементов в моём алгоритме либо что -то изменить в этом чтоб считал точно! Я жду от вас предложений. Заранее спасибо
Тогда в этом случае вообще будет отпадать необходимость в использовании такого медлительного предобуславливателя!Зачем он нужен если он замедляет а не ускоряет!
