Cómo obtener matrices complejas dispersas de mi código a PETSc de manera eficiente

¿Cuál es la forma más eficiente de obtener una matriz dispersa compleja de mi código Fortran a PETSc? Entiendo que esto depende del problema, por lo que traté de dar tantos detalles relevantes como sea posible a continuación.

He estado jugando con el solucionador de valores propios FEAST [1] para problemas del tipo , la dimensión de las matrices y es , y casi todo el tiempo lo paso resolviendo complejo Sistema lineal con M0 a la derecha. N es grande (número de funciones básicas de FE, en 3D), M0 es pequeño (en mi caso, estoy interesado en M0 ~ 20). Las matrices y son reales, simétricas y escasas, y el problema complejo que debe resolverse es , donde$A x = \lambda B x$$A$$B$$N$$N \times N$$A$$B$$zA-B$$z$Es un número complejo. El autor de FEAST parece sugerir que la precisión de la solución a este sistema lineal no tiene que ser muy alta para obtener valores y vectores propios de alta precisión, por lo que algunos solucionadores iterativos rápidos podrían ser una gran solución para esto.

Hasta ahora solo he estado usando Lapack para el sistema complejo, y eso funciona muy bien para en mi computadora. Para más grande , todavía no sé cuál es el mejor solucionador, así que solo quería usar PETSc y jugar con los solucionadores iterativos allí.$N < 1500$$N$

Escribí un controlador C simple, y lo llamo desde Fortran, vea [2] para todo el código, pero el problema es solo con esta parte (actualización: he puesto aquí todas las líneas para crear la matriz, como acabo de darme cuenta que esto es 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);

Lo cual es extremadamente lento (es decir, para N ~ 1500, esto toma tal vez 2 segundos, mientras que resolverlo es inmediato), de hecho, la MatSetValueslínea toma casi el 100% de todo el tiempo para el cálculo del valor propio completo ... La matriz A_ es un Matriz 2D que proviene de Fortran. Traté de desactivar el MAT_IGNORE_ZERO_ENTRIESpero no hizo ninguna diferencia. Entonces, creo que el problema simplemente es que incluso para N moderado como 1500, necesito usar un formato de matriz disperso, ¿es correcto?

Así que implementé el formato CSR en mi código Fortran para las matrices y B (o para la z A - B ) y ahora estoy tratando de descubrir cómo dárselo eficientemente a PETSc. Por ahora, solo quiero que algo funcione secuencialmente, que supere a Lapack. ¿Debo usar la función para esto?$A$$B$$zA-B$MatCreateSeqAIJWithArrays

¿Es esta la forma más eficiente de hacer eso? Dado que las matrices y B no cambian, solo el número complejo z cambia en el algoritmo FEAST, y en un cálculo FE, creo que tanto A como B tienen la misma estructura dispersa, uno probablemente puede mejorar las cosas preasignando el estructura dispersa y luego evaluar rápidamente z A x - B x en cada iteración FEAST ( A x , B x son las matrices de valores del formato CSR), puedo hacer esto fácilmente en Fortran, pero tal vez sea lento llamar siempre$A$$B$$z$$A$$B$$zA_x - B_x$$A_x$$B_x$MatCreateSeqAIJWithArrays, Por lo que podría ser más rápido para hacer todo esto en PETSc vez las matrices y B se transfieren a través.$A$$B$

Me gustaría saber si este enfoque de RSE para este problema es correcto o si lo estoy haciendo mal (claramente mi enfoque original con el MatSetValuesno es óptimo). Gracias por cualquier consejo.

[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

Es importante preasignar correctamente. Esta es casi seguramente la razón por la cual su ensamblaje fue lento. Si está comenzando con una representación de matriz densa, simplemente escanee a través de ella una vez contando el número de nonzeros por fila, luego llame MatSeqAIJSetPreallocation(). Vea estas preguntas frecuentes . La opción MAT_IGNORE_ZERO_ENTRIESrealmente está destinada a ser utilizada cuando hay una leve escasez en esos bloques densos, en lugar de construir la matriz completa en una llamada desde una matriz densa. Por esta razón, no realiza automáticamente la preasignación en función de la escasez de ese bloque.

La creación de una matriz intermedia densa no es escalable en memoria, por lo que eventualmente querrá evitarla. MatSetValues()está realmente destinado a ser utilizado para bloques lógicamente densos en una matriz dispersa. Por lo general, llama una vez por fila (o grupo de filas, típico de los métodos FD) o una vez por elemento (típico de los métodos FEM). Si está traduciendo una matriz dispersa ensamblada existente, solo llame MatSetValues()una vez por fila. Si prefiere omitir la matriz intermedia (mejor rendimiento y menor memoria), solo llame MatSetValues()una vez por elemento.

Tenga en cuenta que puede llamar a PETSc directamente desde Fortran, aunque hay usuarios de todos los dialectos de Fortran entre Fortran 77 básico y las versiones más recientes. Las interfaces se verán bastante crudas para usted como un entusiasta adoptante de las últimas características. Se agradecerán sugerencias sobre formas mantenibles de mejorar el soporte para los últimos dialectos.