Pertanyaan: Metode apa yang tersedia untuk secara akurat dan efisien menghitung struktur sparsity dari matriks elemen hingga?
Info: Saya sedang mengerjakan pemecah Persamaan Tekanan Poisson, menggunakan metode Galerkin dengan basis Lagrange kuadratik, ditulis dalam C, dan menggunakan PETSc untuk penyimpanan matriks jarang dan rutinitas KSP. Untuk menggunakan PETSc secara efisien, saya perlu mengalokasikan memori untuk matriks kekakuan global.
Saat ini, saya sedang melakukan perakitan tiruan untuk memperkirakan jumlah nonzeros per baris sebagai berikut (pseudocode)
int nnz[global_dim]
for E=1 to NUM_ELTS
for i=1 to 6
gi = global index of i
if node gi is free
for j=1 to 6
gj = global index of j
if node gj is free
nnz[i]++
Namun, ini terlalu tinggi nnz karena beberapa interaksi node-node dapat terjadi dalam beberapa elemen.
Saya telah mempertimbangkan untuk mencoba melacak interaksi yang saya temukan, tetapi saya tidak yakin bagaimana melakukan ini tanpa menggunakan banyak memori. Saya juga bisa mengulang simpul, dan menemukan dukungan fungsi dasar yang berpusat pada simpul itu, tetapi kemudian saya harus mencari semua elemen untuk setiap simpul, yang tampaknya tidak efisien.
Saya menemukan pertanyaan baru-baru ini , yang berisi beberapa informasi berguna, terutama dari Stefano M, yang menulis
saran saya adalah untuk mengimplementasikannya dalam python atau C, menerapkan beberapa konsep teori grafik, yaitu mempertimbangkan elemen dalam matriks sebagai tepi dalam grafik dan menghitung struktur sparsity dari matriks adjacency. Daftar daftar atau kamus kunci adalah pilihan umum.
Saya mencari lebih banyak detail dan sumber daya tentang ini. Saya diakui tidak tahu banyak teori grafik, dan saya tidak terbiasa dengan semua trik CS yang mungkin berguna (saya mendekati ini dari sisi matematika).
Terima kasih!
sumber
Jika Anda menetapkan jala Anda sebagai DMPlex dan tata letak data Anda sebagai PetscSection, maka DMCreateMatrix () akan memberi Anda matriks yang telah dialokasikan sebelumnya secara otomatis. Berikut adalah contoh-contoh PETSc untuk Masalah Poisson dan Masalah Stokes .
sumber
Terpilih
Saya pribadi tidak tahu ada cara murah untuk melakukan ini jadi saya hanya melebih-lebihkan angka yaitu, gunakan nilai yang cukup besar untuk semua baris.
Misalnya, untuk mesh yang terstruktur sempurna yang terbuat dari elemen hex 8 node linier, nnzs per baris di blok diagonal dan diagonal off dof * 27 Untuk sebagian besar jerat heks yang sepenuhnya tidak terstruktur secara otomatis, jumlahnya jarang melebihi * 54. Untuk linear tets saya tidak pernah memiliki kebutuhan untuk melampaui dof * 30. Untuk beberapa jerat dengan elemen rasio aspek berbentuk sangat rendah / rendah, Anda mungkin harus menggunakan nilai yang sedikit lebih besar.
Hukumannya adalah konsumsi memori lokal (pada peringkat) adalah antara 2x-5x sehingga Anda mungkin harus menggunakan lebih banyak node komputasi pada cluster Anda daripada biasanya.
Btw saya memang mencoba menggunakan daftar dicari tetapi waktu yang dibutuhkan untuk menentukan struktur sparsity lebih dari perakitan / pemecahan. Tetapi implementasi saya sangat sederhana dan tidak menggunakan informasi tentang edge.
Pilihan lainnya adalah menggunakan rutinitas seperti DMMeshCreateExodus seperti yang ditunjukkan dalam contoh ini .
sumber
Anda ingin menghitung semua koneksi unik (gi, gj), yang menyarankan penempatan semuanya ke dalam wadah asosiatif (non-duplikasi) dan kemudian menghitung kardinalitasnya - dalam C ++ ini akan menjadi std :: set <std :: pair <int, int>>. Di kodesemu, Anda akan mengganti "nnz [i] ++" dengan "s.insert [pair (gi, gj)]", dan kemudian jumlah akhir dari nonzeros adalah s.size (). Ini harus dijalankan dalam waktu O (n-log-n), di mana n adalah jumlah nonzeros.
Karena Anda mungkin sudah mengetahui kisaran gi yang mungkin, Anda dapat "merentang" tabel dengan indeks gi untuk meningkatkan kinerja. Ini mengganti set Anda dengan std :: vector <std :: set <int>>. Anda mengisi itu dengan "v [gi] .insert (gj)", maka jumlah total nonzeros berasal dari menjumlahkan v [gi] .size () untuk semua gi. Ini harus dijalankan dalam waktu O (n-log-k), di mana k adalah jumlah tidak diketahui per elemen (enam untuk Anda - pada dasarnya konstanta untuk sebagian besar kode pde, kecuali jika Anda berbicara tentang metode-hp).
(Catatan - ingin ini menjadi komentar pada jawaban yang dipilih, tetapi terlalu lama - maaf!)
sumber
sumber