Fortran Coder

标题: 关于动态数组的问题 [打印本页]

作者: 不可压缩 时间: 2015-12-5 01:03
标题: 关于动态数组的问题
请教大家一个问题，我有一个很大的稀疏矩阵A(n,n)，我要把所有的非零元素提取出来存入vec，个数未知。
因为尺寸比较大，我想对vec使用动态数组，目前我的方法是。

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

ncnt = 0   !计数器，非零元素的个数

do i = 1,n

 do j = 1,n

    if (A(i,j) /=0)  then

        ncnt =ncnt +1

   endif

 enddo

enddo

allocate(vec(ncnt))



ncnt = 0

do i = 1,n

 do j = 1,n

    if (A(i,j) /=0)  then

        ncnt =ncnt +1

        vec(ncnt) = A(i,j)

   endif

 enddo

enddo

我把A读了两遍，一边用来确定vec的大小，一遍用来赋值。A很大的时候很麻烦，，，
大神们会怎么做呢？

作者: fcode 时间: 2015-12-5 08:44
稀疏矩阵最好使用专门的存储结构和求解器。

如果你真的要自己存储，可以试试这样：

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

program www_fcode_cn

  implicit none

  integer , allocatable :: a(:) , b(:)

  allocate( a(9) )

  a = [1,2,0,0,4,0,0,0,8,0]

  ! 如果你的编译器不支持，请加入以下语句

  !allocate( b( count(a/=0) ) )

  b = pack(a,a/=0)

  write(*,*) b  

end program www_fcode_cn

作者: pasuka 时间: 2015-12-5 09:17
http://people.sc.fsu.edu/~jburkardt/c_src/csparse/csparse.html
Tim Davis写的基于C的稀疏矩阵计算库，也是MATLAB的稀疏矩阵库的早期版本
配合Direct Methods for Sparse Linear Systems, Timothy A. Davis, SIAM, Philadelphia, Sept. 2006.
愿意花半年把这本书啃掉，并拿下C与F混合编程，lz在中国大陆地区的稀疏矩阵编程与计算个人能力排名挤进前25%毫无压力

作者: 不可压缩 时间: 2015-12-7 20:03
我用的就是稀疏矩阵求解器，pardiso。
我要先把稀疏矩阵 A 计算(组装)出来，然后存储成pardiso需要的形式，就是上面说的非零向量vec。

@fcode
什么样的编译器不需要
allocate( b( count(a/=0) ) )
就直接
b = pack(a,a/=0) ?

a 可以是二维吗？

作者: 不可压缩 时间: 2015-12-7 20:05

pasuka 发表于 2015-12-5 09:17
http://people.sc.fsu.edu/~jburkardt/c_src/csparse/csparse.html
Tim Davis写的基于C的稀疏矩阵计算库， ...

没这么大野心啊，我是搞有限元的，求解只是一个工具。
感谢提供的信息

作者: fcode 时间: 2015-12-8 08:30

不可压缩发表于 2015-12-7 20:03
我用的就是稀疏矩阵求解器，pardiso。
我要先把稀疏矩阵 A 计算(组装)出来，然后存储成pardiso需要的形式， ...

我并没有安装所有的编译器，所以也不知道哪些需要，哪些不需要。
但这个用法是 Fortran2008 语法，据我所知，并不是所有编译器都支持到。
你应该自己去尝试，如果无法运行，就加上后面的分配语句。

a 可以是任意维度。

作者: 不可压缩 时间: 2015-12-8 17:51

fcode 发表于 2015-12-8 08:30
我并没有安装所有的编译器，所以也不知道哪些需要，哪些不需要。
但这个用法是 Fortran2008 语法，据我所 ...

哦，我以为是特定的编译器支持这样用。
非常感谢你的回答

作者: pasuka 时间: 2015-12-10 07:24

不可压缩发表于 2015-12-7 20:05
没这么大野心啊，我是搞有限元的，求解只是一个工具。
感谢提供的信息

可惜了，搞不好总刚还是一维变带宽存储或者二维等带宽存储？
CSR、CSC和COO的说明与转换程序，ivf的帮助文档其实写得很明白，主流也就这三种

作者: 不可压缩 时间: 2015-12-15 21:29

pasuka 发表于 2015-12-10 07:24
可惜了，搞不好总刚还是一维变带宽存储或者二维等带宽存储？
CSR、CSC和COO的说明与转换程序，ivf的帮助 ...

有兴趣我们可有交流一下。

我现在的方法是一行一行地组装总刚，依次将该行的非零元素存入上文提到的 vec 中，目前看效率还不错。也有并行的可能。

遗憾的是要提前猜一下总刚里面非零元素有多少个，可能猜大了也可能猜小了

你说的CSR CSC 都是指什么？

作者: pasuka 时间: 2015-12-16 09:02
本帖最后由 pasuka 于 2015-12-16 10:12 编辑

不可压缩发表于 2015-12-15 21:29
有兴趣我们可有交流一下。

我现在的方法是一行一行地组装总刚，依次将该行的非零元素存入上文提到的 vec ...

1、总刚的非零元素个数为啥不能提前算出来？
一维变带宽已经把外圈包络线算出来，包络线内部的零块当然也能算。Bathe和辛克维奇早年的计算机内存太小，再者是LDLT分解不配合AMD这类重排算法的话，零元位置会出现大量非零元，反倒是一维变带宽配合LDLT分解+RCM重排来得实在，但是现在已经扬弃这类办法了。
至于新方法，中文的详见北大的袁明武、陈璞老师以及新科院士田红旗老师的相关论文。此外，数据结构的基础也是必不可缺的，至少得明白各种排序、链表以及hash表
2、CSR、CSC和COO等等，详见ivf的MKL帮助文档的稀疏矩求解章节
btw，多上sourceforge、github以及MATLAB的file exchange上面看看别人写的有限元代码，收获会很多的

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)