2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Неполное разложение Холецкого разреженных матриц
Сообщение13.03.2007, 00:11 


26/11/06
76
В данный момент занимаюсь написанием предобуславливателя к итерационоому решателю 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(неправ) и.т.д. )
Тут надо либо изменить поиск элементов в моём алгоритме либо что -то изменить в этом чтоб считал точно! Я жду от вас предложений. Заранее спасибо

 Профиль  
                  
 
 
Сообщение13.03.2007, 12:41 


13/09/05
153
Москва
Ну так первое, что бросается в глаза - у Вас в IC_Decomposition2 вообще не фигурирует главная диагональ, в то время как в ITL - с элемента главной диагонали начинается каждая строка матрицы.

Во-вторых, можно попробовать оптимизировать Ваш код - зачем искать a[k][j], если a[k][i] == 0 -
Код:
for(k=0;k<i;k++)
{
      m = FindElement(k,i,u,IA,JA);  // ищем элемент a[k][i]
      if (m && r = FindElement(k,JA[j],u,IA,JA))
            sum2 += m*r; //если оба не нулевые то
}

Далее - в двух методах FindElement идет поиск по одной строке, следовательно, можно его объединить. Но ето уже по поводу оптимизации кода в целом, а в Вашем случае это и вовсе не требуется:)

Самое-то главное - зачем искать ненулевые элементы, когда мы специально храним матрицу в удобном для НАС формате, то есть только ненулевые элементы, идущие друг за другом. Следовательно, эту часть неполного преобразования можно свести к простому циклу, без переборов всех возможных вариантов фукцией FindElement, что собственно и сделано в ITL -

Код:
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];
         
}

 Профиль  
                  
 
 
Сообщение13.03.2007, 18:46 


26/11/06
76
Цитата:
Ну так первое, что бросается в глаза - у Вас в IC_Decomposition2 вообще не фигурирует главная диагональ, в то время как в ITL - с элемента главной диагонали начинается каждая строка матрицы.


Когда я сказал , что переписал алгоритм из библиотеки ITL , я забыл упомянуть также то, что я переписал и хранение матриц.
Хранение я переписал так чтобы оно было совместимо с хранением матриц в ITL, т.е. теперь матрица храниться не в 4-х массивах а в 3-х(AN,JA,IA). Теперь все элементы главной диагонали храняться не как раньше в отдельном массиве AD а записываются в AN. Поэтому у меня в методе IC_Decomposition2(double[] AN,int[] JA,int[] IA ) главная диагональ и не фигурирует как отдельная. Всё вроде бы правильно и хранение такое же как в ITL и алгоритм из ITL я переписал один в один, но почему то решение получается не точное!
Ещё раз приведу пример:
Код:
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
//Формируем разреженный формат
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
//Применяем алгоритм, переписанный из ITL
double [] u =  IC_Decomposition2( AN, JA,IA );
//  решение, полученное в рез-те применения  IC_Decomposition2( AN, JA,IA );


u =
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
//Точное решение
u =
3.0
1.0
0.3333333333333333
0.3333333333333333
3.3166247903554
0.6030226891555273
0.30151134457776363
0.6030226891555273
3.1042492870843406
0.5857074126574228
2.750643149492325
0.24236755930688816
3.439498052781032
0.2584356424246762
2.8284271247461903
2.731018774006702     


Цитата:
Во-вторых, можно попробовать оптимизировать Ваш код - зачем искать a[k][j], если a[k][i] == 0


Да. И в правду зачем.....
Код:
for(k=0;k<i;k++) {
                  m=FindElement(k,i,u,IA,JA);  // ищем элемент a[k][i]
              if (m!=0){
                  r=FindElement(k,JA[j],u,IA,JA); // ищем элемент a[k][j]   
                 if (r!=0)  sum2+=m*r;//если оба не нулевые то перемножаем
              }                                     
          }   


Несмотря на эту небольшую оптимизацию, мною написанный алгоритм всё ещё очень медлителен.
Цитата:
Самое-то главное - зачем искать ненулевые элементы, когда мы специально храним матрицу в удобном для НАС формате, то есть только ненулевые элементы, идущие друг за другом.

А как мне ещё по - другому узнать есть ли, ненулевой элемент в i-ом столбце k-ой строки как не перебором всех ненулевых элементов k-ой строки?

 Профиль  
                  
 
 А должны ли разложения быть одинаковыми, если неполные?
Сообщение13.03.2007, 23:39 


03/09/05
217
Bulgaria
Гипотеза?
В этом Форуме приводилась ссылка на учебное пособие "Методы решения СЛАУ большой размерности", авторы М. Ю. Баландин и Э. П. Шурина.
Там на стр. 12-17 отлично показано, что неполное LU разложение (ILU), потому что неполное, и вовсе не должно равняться по всему множеству пар индексов исходной матрицы, ее элементам.
Но зато - масштабирующие множители легко вычислимы и легко обратимы.
И масштабированная матрица точнее обрабатывается далее, а потом расмаштабируется.
И наверное, два разные способа неполного разложения, вполне возможно будут приводить к разным, хотя и близким ввиду близости к исходной матрицы, множителям разложений.

 Профиль  
                  
 
 
Сообщение14.03.2007, 14:43 


13/09/05
153
Москва
Цитата:
Несмотря на эту небольшую оптимизацию, мною написанный алгоритм всё ещё очень медлителен.


Не знаю, это как посмотреть:)
Воспроизвел Ваш код и меня получилось следующее:
1. Исходный код при следующем вызове
Код:
for (int i1 = 0; i1 < 100000; i1++)
      IC_Decomposition(nANSize, nRows, AN, AD, JA, IA);

работает в среднем 0,2018 сек.

2. Немного модифицированный, где
Код:
m = FindElement(k, i, u, IA, JA);  // ищем элемент a[k][i]
            if (m != 0 && (r = FindElement(k, JA[j], u, IA, JA)) != 0)
               sum2 += m*r;

работает 0,17617 сек.

3. Переделал поиск суммы, таким образом функция IC_Decomposition имеет вид
Код:
   for (i = 1; i < n; i++)
   {
      udiag[i] = sqrt(AD[i] - sum[i]);//находим диагональные элементы
      for (int j = IA[i]; j < IA[i + 1]; j++)
      {
         u[j] = AN[j];
         for (int k = 0; k < i; k++)
         {
            int n1 = -1;
            for (int ii = IA[k]; 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] /= udiag[i];
         sum[JA[j]] += u[j]*u[j];
      }
   }

Этот код работает примерно 0,0996 сек.

Теперь про ITL - у меня он работает примерно так же - 0,1067 сек и это в прицнипе правильно:), так как у них, и в моей версии функции есть 4 цикла и несколько проверок без лишних вызовов функций и лишних действий.

Цитата:
А как мне ещё по - другому узнать есть ли, ненулевой элемент в i-ом столбце k-ой строки как не перебором всех ненулевых элементов k-ой строки?

Я имел ввиду, что элементы одной строки идут подряд, Вы ищите a[k][i] и a[k][j], при этом i < j, таким образом поиск можно объединить. Если мы уже рассматриваем j-элемент, то прекращаем поиск. Если попутно нашли i-элемент - вычисляем сумму. Не нужно искать сами ненулевые элементы, нужно найти их сумму:)

 Профиль  
                  
 
 
Сообщение15.03.2007, 13:42 


26/11/06
76
Цитата:
Воспроизвел Ваш код и меня получилось следующее:
1. Исходный код при следующем вызове
Код:
Код:
for (int i1 = 0; i1 < 100000; i1++)
      IC_Decomposition(nANSize, nRows, AN, AD, JA, IA);


работает в среднем 0,2018 сек.



А в качестве массивов AN,AD,JA,IA вы какие брали и какой размерности для теста? С моего 1-ого примера:
Код:
Разреженная матрица ------------
Формат хранения: 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
-----------------------------------------------

или свои?


Цитата:
Теперь про ITL - у меня он работает примерно так же - 0,1067 сек и это в прицнипе правильно, так как у них, и в моей версии функции есть 4 цикла и несколько проверок без лишних вызовов функций и лишних действий.

Так почему же всё таки мною переписанный код из ITL даёт не точное решение?

 Профиль  
                  
 
 
Сообщение15.03.2007, 14:12 


13/09/05
153
Москва
Цитата:
А в качестве массивов AN,AD,JA,IA вы какие брали и какой размерности для теста? С моего 1-ого примера:

Конечно, Ваши данные:)

Цитата:
Так почему же всё таки мною переписанный код из ITL даёт не точное решение?

А Вы посмотрите внимательно, там немного по другому считается:)
Бросается в глаза , что у Вас диагональный элемент udiag[i] = sqrt(AD[i] - sum[i]);
а в ITL - AN_[d] = sqrt(AN_[d]);
Хотя тут нужно проверить - они ведь исходную матрицу AN подправляют на месте.

При тестировании я убрал выделение памяти из кода IC_Decomposition, так как IC_Decomposition2 не выделяет памяти, а пишет поверх матрицы AN.
В моей реализации IC_Decomposition пишет в глобальные массивы u, udiag. sum.
Таким образом получилось сравнить работу двух алгоритмов.

В Вашем исходном коде есть выделение памяти в теле функции и его запускать в цикле смысла нету, большая часть времени тратится на выделение памяти:)

 Профиль  
                  
 
 
Сообщение17.03.2007, 20:16 


26/11/06
76
При увеличении размерности решаемой задачи мой алгоритм начинает существенно проигрывать по скорости ITL - вскому.

Тесты (на JDK 5.0):

Матрица 7x7 (та которую я уже ранее приводил в качестве примера):
При таком вызове:
Код:
for (int i1 = 0; i1 < 100000; i1++) IC_Decomposition( AN, AD, JA, IA);

мой алгоритм, с переделанным вами поиском суммы , выполняется
0,063 сек.
Алгоритм IC_Decomposition2, переписанный мною из ITL выполняется
за 0,047 сек.
Вроде бы почти одинаково. Но стоит увеличить размерность и.....
Я взял в качестве след. теста матрицу размером 30x30.

Код:
46 -11 0 -21 8 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-11 6 -2 8 -4 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 -2 10 -9 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-21 8 -9 37 -19 -1 -16 11 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
8 -4 3 -19 14 -2 11 -10 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
9 -3 2 -1 -2 10 -7 5 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -16 11 -7 26 -22 -1 -10 11 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 11 -10 5 -22 26 -1 11 -16 -7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 7 -5 2 -1 -1 10 -5 7 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 -10 11 -5 14 -19 -2 -4 8 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 11 -16 7 -19 37 -1 8 -21 -9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 5 -7 2 -2 -1 10 -3 9 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 -4 8 -3 6 -11 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 8 -21 9 -11 46 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 3 -9 2 -2 0 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 43 -106 22 -36 78 33 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -106 413 5 78 -192 -81 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 22 5 91 -33 81 22 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -36 78 -33 121 -184 19 -85 106 53 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 78 -192 81 -184 335 11 106 -142 -69 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 33 -81 22 19 11 91 -53 69 22 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -85 106 -53 228 -213 16 -142 106 69 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 106 -142 69 -213 228 16 106 -85 -53 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 53 -69 22 16 16 91 -69 53 22 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -142 106 -69 335 -184 11 -192 78 81
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 106 -85 53 -184 121 19 78 -36 -33
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 69 -53 22 11 19 91 -81 33 22
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -192 78 -81 413 -106 5
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 78 -36 33 -106 43 22
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 81 -33 22 5 22 91

Упаковал её в 4 - массива, вот они:
Код:
------------- Разреженная матрица ------------
Формат хранения: CSIR
AN:-11.0 -21.0 8.0 9.0 -2.0 8.0 -4.0 -3.0 -9.0 3.0 2.0 -19.0 -1.0 -16.0 11.0 7.0 -2.0 11.0 -10.0 -5.0 -7.0 5.0 2.0 -22.0 -1.0 -10.0 11.0 5.0 -1.0 11.0 -16.0 -7.0 -5.0 7.0 2.0 -19.0 -2.0 -4.0 8.0 3.0 -1.0 8.0 -21.0 -9.0 -3.0 9.0 2.0 -11.0 -2.0 -106.0 22.0 -36.0 78.0 33.0 5.0 78.0 -192.0 -81.0 -33.0 81.0 22.0 -184.0 19.0 -85.0 106.0 53.0 11.0 106.0 -142.0 -69.0 -53.0 69.0 22.0 -213.0 16.0 -142.0 106.0 69.0 16.0 106.0 -85.0 -53.0 -69.0 53.0 22.0 -184.0 11.0 -192.0 78.0 81.0 19.0 78.0 -36.0 -33.0 -81.0 33.0 22.0 -106.0 5.0 22.0
JA:1 3 4 5 2 3 4 5 3 4 5 4 5 6 7 8 5 6 7 8 6 7 8 7 8 9 10 11 8 9 10 11 9 10 11 10 11 12 13 14 11 12 13 14 12 13 14 13 14 16 17 18 19 20 17 18 19 20 18 19 20 19 20 21 22 23 20 21 22 23 21 22 23 22 23 24 25 26 23 24 25 26 24 25 26 25 26 27 28 29 26 27 28 29 27 28 29 28 29 29
IA:0 4 8 11 16 20 23 28 32 35 40 44 47 49 49 49 54 58 61 66 70 73 78 82 85 90 94 97 99 100 100
AD: 46.0 6.0 10.0 37.0 14.0 10.0 26.0 26.0 10.0 14.0 37.0 10.0 6.0 46.0 10.0 43.0 413.0 91.0 121.0 335.0 91.0 228.0 228.0 91.0 335.0 121.0 91.0 413.0 43.0 91.0
-----------------------------------------------

Решил её также 100.000 раз модифицированным вами IC_Decomposition ( AN, AD, JA, IA);
Время решения: 5.187 сек :(

Упаковал её в 3-массива специально для решения под IC_Decomposition2 (AN,JA,IA):
Код:
------------- Разреженная матрица ------------
Формат хранения: CSIR
AN: 46.0 -11.0 -21.0 8.0 9.0 6.0 -2.0 8.0 -4.0 -3.0 10.0 -9.0 3.0 2.0 37.0 -19.0 -1.0 -16.0 11.0 7.0 14.0 -2.0 11.0 -10.0 -5.0 10.0 -7.0 5.0 2.0 26.0 -22.0 -1.0 -10.0 11.0 5.0 26.0 -1.0 11.0 -16.0 -7.0 10.0 -5.0 7.0 2.0 14.0 -19.0 -2.0 -4.0 8.0 3.0 37.0 -1.0 8.0 -21.0 -9.0 10.0 -3.0 9.0 2.0 6.0 -11.0 -2.0 46.0 10.0 43.0 -106.0 22.0 -36.0 78.0 33.0 413.0 5.0 78.0 -192.0 -81.0 91.0 -33.0 81.0 22.0 121.0 -184.0 19.0 -85.0 106.0 53.0 335.0 11.0 106.0 -142.0 -69.0 91.0 -53.0 69.0 22.0 228.0 -213.0 16.0 -142.0 106.0 69.0 228.0 16.0 106.0 -85.0 -53.0 91.0 -69.0 53.0 22.0 335.0 -184.0 11.0 -192.0 78.0 81.0 121.0 19.0 78.0 -36.0 -33.0 91.0 -81.0 33.0 22.0 413.0 -106.0 5.0 43.0 22.0 91.0
JA: 0 1 3 4 5 1 2 3 4 5 2 3 4 5 3 4 5 6 7 8 4 5 6 7 8 5 6 7 8 6 7 8 9 10 11 7 8 9 10 11 8 9 10 11 9 10 11 12 13 14 10 11 12 13 14 11 12 13 14 12 13 14 13 14 15 16 17 18 19 20 16 17 18 19 20 17 18 19 20 18 19 20 21 22 23 19 20 21 22 23 20 21 22 23 21 22 23 24 25 26 22 23 24 25 26 23 24 25 26 24 25 26 27 28 29 25 26 27 28 29 26 27 28 29 27 28 29 28 29 29
IA: 0 5 10 14 20 25 29 35 40 44 50 55 59 62 63 64 70 75 79 85 90 94 100 105 109 115 120 124 127 129 130
-----------------------------------------------


Решил её 100.000 раз IC_Decomposition2 (AN,JA,IA)
Время решения составило всего 0.828 сек.
Как вы видите IC_Decomposition2 практически в 6 раз быстрее решает чем IC_Decomposition.
А при ещё большем увеличении размерности задачи этот разрыв будет ещё больше увеличиваться. А мне в дальнейшем потребуется решать задачи с матрицей 100.000x100.000 и более и тогда, если использовать IC_Decomposition - предобуславливатель, очень много времени будет уходить только на предобуславливание :( Тогда в этом случае вообще будет отпадать необходимость в использовании такого медлительного предобуславливателя!Зачем он нужен если он замедляет а не ускоряет!


Цитата:

А Вы посмотрите внимательно, там немного по другому считается
Бросается в глаза , что у Вас диагональный элемент udiag[i] = sqrt(AD[i] - sum[i]);
а в ITL - AN_[d] = sqrt(AN_[d]);
Хотя тут нужно проверить - они ведь исходную матрицу AN подправляют на месте.



Диагональный элемент udiag[i] = sqrt(AD[i] - sum[i]) у меня считается так только в IC_Decomposition, придуманным мною. А в IC_Decomposition2 диагональный элемент считается точно так же как и у них (в ITL) и хранение у меня организовано точно так же как и в ITL (3 -массива).(Посмотрите ещё раз код IC_Decomposition2 и код из ITL в 1-ом посте,принимая во в нимание то что val = = AN , indx = = JA и pntr = = IA.Один в один переписал!!!). К сожалению IC_Decomposition2 быстр но неточен :cry:

 Профиль  
                  
 
 
Сообщение19.03.2007, 13:37 


13/09/05
153
Москва
Проверил - действительно ITL шустрее - у меня это 6.97 против 1.26 секунд для матрицы 30х30.

Цитата:
Диагональный элемент udiag[i] = sqrt(AD[i] - sum[i]) у меня считается так только в IC_Decomposition, придуманным мною. А в IC_Decomposition2 диагональный элемент считается точно так же как и у них (в ITL) и хранение у меня организовано точно так же как и в ITL (3 -массива).(Посмотрите ещё раз код IC_Decomposition2 и код из ITL в 1-ом посте,принимая во в нимание то что val = = AN , indx = = JA и pntr = = IA.Один в один переписал!!!).


Здесь я имел ввиду, что в Вашем варианте (IC_Decomposition) и в ITL разложение считается немного по-разному. Просто выпишите алгоритм из ITL и сравните то, что делают они в ITL, с тем, что есть в Вашей книге (источнике).

Далее - они не используют новых массивов - а пишут поверх. Это избавляет от присвоений.
Диагональные элементы считаются по-другому - просто как корень, а не как
udiag[i] = sqrt(AD[i] - sum[i]) - следовательно, нет одной операции вычитания, а также не используется значение суммы (и сам массив sum), следовательно не вычисляется сумма и не делается присвоение вида sum[JA[j]] += u[j]*u[j] с произведением.
В остальном все практически также - 4 цикла, одна операция корня и одно умножение и одно деление.

Следует отметить, что главная цель предобуславливателя - ускорить сходимость матрицы и процесс вычисления. Чтоб у Вас не было иллюзий - матрица 100,000х100,000 - считатается несколько минут, поэтому если в ITL разложение делатся за 1 сек, а у Вас - за 5 сек, то это не страшно.

Цитата:
К сожалению IC_Decomposition2 быстр но неточен
- у них просто другой алгоритм разложения, Вы по решению проверьте. Я пользуюсь ITL уже давно и пока не жалуюсь. Само по себе преобразование никому не нужно, важно решение:))

Советую протестировать на решении конкретных задач. Здесь я запускал в цикле матрицу малых размеров, а по хорошему нужно взять реальную МКЭ-задачу с большим числом конечных элементов, порядка 100,000 и более. 100,000 КЭ первого порядка - это примерно 50,000 узлов, для КЭ второго порядка - 200,000.
(В общем случае, число узлов КЭ первого порядка равно половине числа элементнов, число узлов для КЭ второго порядка - 2х число элементов).
Расчитать, посмотреть что-да-как:)
При этом можно посмотреть на число итераций решателя (Вы же итерационный собираетесь использовать, CG) для Вашего преобразования и их преобразования. Может быть после Вашего преобразования быстрее считаться будет:)))

 Профиль  
                  
 
 
Сообщение29.03.2007, 16:46 


26/11/06
76
Цитата:
Здесь я имел ввиду, что в Вашем варианте (IC_Decomposition) и в ITL разложение считается немного по-разному. Просто выпишите алгоритм из ITL и сравните то, что делают они в ITL, с тем, что есть в Вашей книге (источнике).


Да IC_Decomposition и ITL совсем разное. Но IC_Decomposition2 и ITL один и тот же.
Цитата:
Я пользуюсь ITL уже давно и пока не жалуюсь

А какое решение той самой первой матрицы 7x7 давал ITL-овское разложение у вас, можете выложить?
Цитата:

у них просто другой алгоритм разложения, Вы по решению проверьте.
При этом можно посмотреть на число итераций решателя (Вы же итерационный собираетесь использовать, CG) для Вашего преобразования и их преобразования. Может быть после Вашего преобразования быстрее считаться будет


Т.е вы считаете то что у меня IC_Decomposition2 считает немного не точно это так и должно быть? Что у них (в ITL) после разложения получается тоже такое же не совсем точное решение, приближенное к точному неполному разложению Холецкого.


Протестировал на реальной задачи, содержащей 2000 узлов и 1999 кон. элементов. (Арка). Матрица жесткости получилась размером 6000 x 6000. Число ненулевых элементов составило 89 000.

Протестировал эту задачу на трёх ситуациях:


Без предобуславливателя. CG сошелся за 0,078 сек на 27-ой итерации.
С предобуславливателем IC_Decomposition. Само разложение делалось около 20 секунд.(достаточно долго как вы видите). CG сошелся за 1 - ну итерацию за 0,031 сек.
С itl - овским предобуславливателем IC_Decomposition2. Разложение делалось 0,031 сек. CG сошелся на 9 – ой итерации за 0.062 секунд.

Вывод: Если использовать в качестве предобуславливателя IC_Decomposition (мой) , то уйма времени решения всей задачи уходит только на само разложение, зато этот алгоритм обеспечивает быструю сходимость CG.
Если использовать в качестве предобуславливателя itl - овский IC_Decomposition2, то время на разложение практически не затрачивается , но из – за того что это разложение отличается от точного неполного разложения Холецкого процесс сходимости CG затягивается. Вот если бы сделать так чтобы IC_Decomposition2 давал точное разложение то CG сошелся бы за 1 итерацию и само разложение заняло бы очень мало времени , т.е ускорился бы процесс решения всей задачи.

При тестировании нескольких реальных задач наткнулся на очень негативный эффект Неполного разложения Холецкого. На некоторых плохо обусловленных симметричных положительно определенных матрицах неполное разложение Холецкого (как моё так и Itl - овское) вообще не выполняется.Это происходит из – за того что в мы IC жестко пренебрегаем появлением новых ненулевых элементов в тех позициях матрицы, где раньше стояли нули , то при вычислении диагональных элементов значение подкоренного выражения становится намного меньше того , если бы мы учитывали все новые появившиеся ненулевые элементы. И поэтому во многих задачах под корнем получаются отрицательные значения и такое разложение не выполняется. Что делать в таких ситуациях???

Вот пример такой матрицы
http://rapidshare.com/files/23344051/notSolveMatrix.rar
. В архиве 7 файлов. AN.txt,JA.txt,IA.txt,AD.txt – это эта матрица представленная для IC_Decomposition. AN2.txt,JA2.txt и IA2.txt – это матрица, представленная для itl - овского IC_Decomposition2. Не 1-ый не 2-ой алгоритм её не разложил. Попробуйте у себя её разложить неполным разложением Холецкого.
Сначала я подумал что она не положительно определена. Но после то го как я её конвертировал в полное представление (со всеми нулями) и сделал полное разложение Холецкого то мои сомнения по поводу не положительно определенности этой матрицы отпали. Вектор полученных неизвестных по методу полного Холецкого я записал в 7-ой файл Вектор решений.txt.

Цитата:
Советую протестировать на решении конкретных задач


Можете вы мне сгенерировать в вашей программе несколько матриц и векторов для реальных задач (можно в полном виде).
Размерности: 10000x10000,30000x30000,50000x50000,70000x70000,100000x100000 и 200000x200000.

 Профиль  
                  
 
 
Сообщение31.03.2007, 14:06 


13/09/05
153
Москва
Цитата:
Т.е вы считаете то что у меня IC_Decomposition2 считает немного не точно это так и должно быть? Что у них (в ITL) после разложения получается тоже такое же не совсем точное решение, приближенное к точному неполному разложению Холецкого.

Цель предобуславливателя - ускорить процесс схождения, и в ITL это делается:)
Там кстати есть еще и modified_cholesky - беглый просмотр по коду говорит, что там немного по другому считатется, но нужно детально посмотреть:)

Цитата:
CG сошелся бы за 1 итерацию и само разложение заняло бы очень мало времени , т.е ускорился бы процесс решения всей задачи.

За 1 итерацию это значит, что всю работу делает разложение, а не решатель.
На то он и итерационный, чтоб решать за несколько итераций, при этом ему можно немного помочь, применив предобуславливание.

Цитата:
При тестировании нескольких реальных задач наткнулся на очень негативный эффект Неполного разложения Холецкого.

У меня для ITL за долгое время его использования такого замечено не было.
Попробовал открыть Вашу матрицу (AN2, IA2, JA2) - там у Вас напутано с IA2, JA2.
Там есть нули в позициях где их быть не должно:)
Поэтому проверить не получилось, попробуйте занова создать их и выложить.

Цитата:
Можете вы мне сгенерировать в вашей программе несколько матриц и векторов для реальных задач (можно в полном виде).
Размерности: 10000x10000,30000x30000,50000x50000,70000x70000,100000x100000 и 200000x200000.

Можно, но ето целая история.
Самый простой вариант - это приделать Вам генератор сетки и некоторое средство для задания геометрии и исходных данных. Такие задачи нужно делать не по частям, а в купе -
Задание исходных данных-Генератор сетки-Создание СЛАУ-Решатель-Построение графиков.
А потом все отлаживать каждую часть, что намного проще чем делать это в слепую.

Из доступных генератор сетки могу посоветовать GeomPack++, имеет достаточно простой формат входных и выходных данных, таким образом нужно всего лишь сделать генерации входных файлов, запуск программы GeomPack и чтение сетки КЭ.
Сетка в GeomPack++ для 2D кривая (и Делоне, и Advancing Front), но для отладки вполне подойдет. В нем кстати можно задавать число суммарное конечных элементов в сетке.

 Профиль  
                  
 
 
Сообщение31.03.2007, 16:41 


26/11/06
76
Цитата:
Попробовал открыть Вашу матрицу (AN2, IA2, JA2) - там у Вас напутано с IA2, JA2.
Там есть нули в позициях где их быть не должно
Поэтому проверить не получилось, попробуйте занова создать их и выложить.


Да действительно JA2 был неправильно запсан.Вот новые файлы.
http://rs33.rapidshare.com/files/23661269/4.rar

 Профиль  
                  
 
 
Сообщение02.04.2007, 13:56 


13/09/05
153
Москва
Попробовал разные решатели и предобуславливатели - процесс не сходится.
Матрица плохо обусловлена.
Возник такой - Вы ничего в создании матрицы не напутали?
Судя по матрице - неизвестно ни одно ГУ первого рода. Вектор правых частей весь одинаковый, что тоже не понятно.
Просто у меня при решении уравнения Пуассона такой ерунды с решателем никогда не возникало, все решается нормально:)

 Профиль  
                  
 
 
Сообщение02.04.2007, 15:54 


26/11/06
76
Цитата:
Попробовал разные решатели и предобуславливатели - процесс не сходится.

У меня сходится после 73659 итераций по методу CG без использования предобуславливателя за 31,5 сек.
А какие вы пробовали ?
Цитата:
Возник такой - Вы ничего в создании матрицы не напутали?

Нет.Матрица создавалась в МКЭ программе расчета плоских стержневых конструкций. В качестве схемы была выбрана арка состоящая из 1000 узлов и 999 конечных элементов. В каждом узле 3 степени свободы.
Цитата:
Вектор правых частей весь одинаковый, что тоже не понятно

Вектор правых частей можете поставить любой. От него не зависит само разложение.

Вы мне самого главного не сказали, получилось ли у вас разложить данную матрицу Неполным Холецким , взятым из ITL?

 Профиль  
                  
 
 
Сообщение02.04.2007, 20:39 


13/09/05
153
Москва
Цитата:
Вы мне самого главного не сказали, получилось ли у вас разложить данную матрицу Неполным Холецким, взятым из ITL?

Да вроде разложилась:)

cholesky + cg - делает много итераций и не сходится
modified_cholesky + cg - здесь на второй итерации все рушится, вот тут походу и не разложилась:)
ilu +cg - тоже не сходится

Кроме того, проверил по IC_Decomposition2 - нулевых подкоренных значений там нет.

Цитата:
У меня сходится после 73659 итераций по методу CG без использования предобуславливателя за 31,5 сек.

73659 - это уже перебор. Так обычно максимум за 1000 итераций все сходится. У меня даже ограничение в программе стояло на 1000, если не сошлось - то значит и не сойдется:)
Потом с какой точностью сошлось - относительной/абсолютной.
Матрица такой размерности должна решаться за десятые доли секунды.

Как пример - провел расчет эл. поля - сетка 117136 КЭ второго порядка, одна степень свободы в узле, полное число узлов - примерно 2*117136.
При задании относительной точности 1е-8 (а больше и нужно) связка cholesky + cg -
148 iterations
0.512286 is actual final residual.
9.59006e-009 is actual relative tolerance achieved.

Цитата:
Вектор правых частей можете поставить любой. От него не зависит само разложение.

Это да, но я имел ввиду то, что непонятно что Вы решаете.

Еще хороший вариант проверить решение Вашей матрицы - это MATLAB.
В нем прямо из файлов прочитать матрицы и попробовать разные решатели для этой матрицы.

И, кстати, почему Холецкий, когда он вроде Холесский?:)))

Подводя итоги -
1. Программу МКЭ нужно отлаживать в целом, что можно просто сделать самому - сделать, а то что сложно - взять на первое время готовое. Потом при наличии свободного времени и желания сделать свое (генератор сетки, задание геометрии и построение графиков).
2. Когда у Вас будет в том или ином виде препроцессор и постпроцессор - тогда и отлаживать решатель на конкретных задачах.
3. Основная задача предобуславливателя - облегчить работу итерационного решателя, чтоб он немного быстрее сошелся.
4. Процесс решения можно проверить на задачах с известным аналитическим решением. Сам решатель - сравнить с решением готовыми библиотеками (тот же ITL) или с МАТЛАБом.
5. Большая часть коммерческого софта использует связку cholesky+CG.
В FEMLAB - если не ошибаюсь, по умолчанию это прямой решатель.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 24 ]  На страницу 1, 2  След.

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group