Como obter matrizes complexas esparsas do meu código para o PETSc com eficiência

Qual é a maneira mais eficiente de obter uma matriz esparsa complexa do meu código Fortran para o PETSc? Entendo que isso depende do problema, então tentei fornecer o maior número possível de detalhes relevantes abaixo.

Eu estou brincando com o solucionador de autovalor FEAST [1] para problemas do tipo , a dimensão das matrizes e é , e praticamente todo o tempo gasto resolvendo sistema linear complexo com o lado direito M0. N é grande (número de funções básicas da FE, em 3D), M0 é pequeno (no meu caso, estou interessado em M0 ~ 20). As matrizes e são reais, simétricas e esparsas, e o problema complexo que precisa ser resolvido é , onde$A x = \lambda B x$$A$$B$$N$$N \times N$$A$$B$$zA-B$$z$é um número complexo. O autor do FEAST parece sugerir que a precisão da solução para esse sistema linear não precisa ser muito alta para obter valores próprios e vetores próprios precisos, portanto alguns solucionadores iterativos rápidos podem ser uma ótima solução para isso.

Até agora, eu apenas uso o Lapack para o sistema complexo, e isso funciona muito bem para no meu computador. Para maior , ainda não sei qual é o melhor solucionador, então eu queria apenas usar o PETSc e brincar com os solucionadores iterativos lá.$N < 1500$$N$

Eu escrevi um driver C simples, e chamo-o de Fortran, veja [2] para todo o código, mas o problema é apenas com essa parte (atualização: coloquei aqui todas as linhas para criar a matriz, como acabei de perceber que isso é relevante):

ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
for (i=0; i<n; i++) col[i] = i;
ierr = MatSetValues(A,n,col,n,col,A_,INSERT_VALUES);CHKERRQ(ierr);

O que é extremamente lento (ou seja, para N ~ 1500, isso leva talvez 2s, enquanto a resolução é imediata); na verdade, a MatSetValueslinha leva praticamente 100% de todo o tempo para todo o cálculo do autovalor ... A matriz A_ é uma Matriz 2D que vem do Fortran. Eu tentei desativar o MAT_IGNORE_ZERO_ENTRIESmas não fez nenhuma diferença. Então, acho que o problema simplesmente é que, mesmo para N moderado como 1500, preciso usar algum formato de matriz esparsa, está correto?

Portanto, implementei o formato CSR no meu código Fortran para as matrizes e (ou para o ) e agora estou tentando descobrir como fornecê-lo ao PETSc de forma eficiente. Por enquanto, eu só quero que algo funcione sequencialmente, que supera Lapack. Devo usar a função para isso?$A$$B$$zA-B$MatCreateSeqAIJWithArrays

Essa é a maneira mais eficiente de fazer isso? Como as matrizes e não mudam, apenas o número complexo muda no algoritmo FEAST e, em um cálculo de FE, acho que ambos e têm a mesma estrutura esparsa, provavelmente é possível melhorar ainda mais as coisas pré-alocando o estrutura esparsa e, em seguida, apenas avalie rapidamente em cada iteração FEAST (o , são as matrizes de valores do formato CSR), posso fazer isso facilmente no Fortran, mas talvez seja lento sempre chamar , portanto, pode ser mais rápido para fazer tudo isso no PETSc uma vez que as matrizes e$A$$B$$z$$A$$B$$zA_x - B_x$$A_x$$B_x$MatCreateSeqAIJWithArrays$A$$B$ são transferidos.

Gostaria de saber se essa abordagem de RSE para esse problema está correta ou se estou fazendo errado (claramente minha abordagem original com a opção MatSetValuesnão é ótima). Obrigado por todas as dicas.

[1] http://www.ecs.umass.edu/~polizzi/feast/

[2] https://github.com/certik/hfsolver/pull/14.

[3] http://www.mcs.anl.gov/petsc/petsc-3.1/docs/manualpages/Mat/MatCreateSeqAIJWithArrays.html

É importante pré-alocar corretamente. Esta é quase certamente a razão pela qual sua montagem foi lenta. Se você estiver começando com uma representação de matriz densa, basta digitalizá-la uma vez, contando o número de não-zeros por linha, e ligue MatSeqAIJSetPreallocation(). Veja este FAQ . A opção MAT_IGNORE_ZERO_ENTRIESrealmente deve ser usada quando houver uma leve dispersão nesses blocos densos, em vez de construir toda a matriz em uma chamada a partir de uma matriz densa. Por esse motivo, ele não realiza a pré-localização automaticamente com base na escassez desse bloco.

Criar uma matriz intermediária densa não é escalável em memória; portanto, você desejará evitá-la. MatSetValues()é realmente destinado a ser usado para blocos logicamente densos em uma matriz esparsa. Você normalmente chama uma vez por linha (ou grupo de linhas, típico dos métodos FD) ou uma vez por elemento (típico dos métodos FEM). Se você estiver traduzindo uma matriz esparsa montada existente, basta ligar MatSetValues()uma vez por linha. Se você preferir pular a matriz intermediária (melhor desempenho e menos memória), basta chamar MatSetValues()uma vez por elemento.

Observe que você pode ligar para o PETSc diretamente do Fortran, embora existam usuários de todos os dialetos do Fortran entre o Fortran 77 básico e as versões mais recentes. As interfaces parecerão um pouco grosseiras para você como um adotante ansioso dos recursos mais recentes. Seriam apreciadas sugestões de maneiras sustentáveis ​​de melhorar o suporte para os dialetos mais recentes.