不可压缩 发表于 2015-12-5 01:03:27

关于动态数组的问题

请教大家一个问题，我有一个很大的稀疏矩阵A(n,n)， 我要把所有的非零元素提取出来存入vec，个数未知。
因为尺寸比较大，我想对vec使用动态数组，目前我的方法是。
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:47

稀疏矩阵最好使用专门的存储结构和求解器。

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

program www_fcode_cn
implicit none
integer , allocatable :: a(:) , b(:)
allocate( a(9) )
a =
! 如果你的编译器不支持，请加入以下语句
!allocate( b( count(a/=0) ) )
b = pack(a,a/=0)
write(*,*) b
end program www_fcode_cn

pasuka 发表于 2015-12-5 09:17:40

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:27

我用的就是稀疏矩阵求解器，pardiso。
我要先把稀疏矩阵 A 计算(组装)出来，然后存储成pardiso需要的形式，就是上面说的非零向量vec。

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

a 可以是二维吗？

不可压缩 发表于 2015-12-7 20:05:04

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:05

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

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

a 可以是任意维度。

不可压缩 发表于 2015-12-8 17:51:46

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

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

pasuka 发表于 2015-12-10 07:24:49

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

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

不可压缩 发表于 2015-12-15 21:29:09

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

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

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

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

你说的CSR CSC 都是指什么？

pasuka 发表于 2015-12-16 09:02:28

本帖最后由 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上面看看别人写的有限元代码，收获会很多的
页: [1]
查看完整版本: 关于动态数组的问题