petsc4py: Creating AIJ Matrix from csc_matrix resu

问题:

I am trying to create a petsc-matrix form an already existing csc-matrix. With this in mind I created the following example code:

import numpy as np
import scipy.sparse as sp
import math as math
from petsc4py import PETSc
n=100

A = sp.csc_matrix((n,n),dtype=np.complex128)
print A.shape
A[1:5,:]=1+1j*5*math.pi

p1=A.indptr
p2=A.indices
p3=A.data
petsc_mat = PETSc.Mat().createAIJ(size=A.shape,csr=(p1,p2,p3))

This works perfectly well as long as the matrix only consist of real values. When the matrix is complex running this piece of code results in a TypeError: Cannot cast array data from dtype('complex128') to dtype('float64') according to the rule 'safe'. I tried to figure out where the error occurs exactly, but could not make much sense of the traceback:

petsc_mat = PETSc.Mat().createAIJ(size=A.shape,csr=(p1,p2,p3))  File "Mat.pyx", line 265, in petsc4py.PETSc.Mat.createAIJ (src/petsc4py.PETSc.c:98970)
File "petscmat.pxi", line 662, in petsc4py.PETSc.Mat_AllocAIJ (src/petsc4py.PETSc.c:24264)
File "petscmat.pxi", line 633, in petsc4py.PETSc.Mat_AllocAIJ_CSR (src/petsc4py.PETSc.c:23858)
File "arraynpy.pxi", line 136, in petsc4py.PETSc.iarray_s (src/petsc4py.PETSc.c:8048)
File "arraynpy.pxi", line 117, in petsc4py.PETSc.iarray (src/petsc4py.PETSc.c:7771)

Is there an efficient way of creating a petsc matrix (of which i want to retrieve some eigenpairs later) from a complex scipy csc matrix ?

I would be really happy if you guys could help me find my (hopefully not too obvious) mistake.

回答1:

I had troubles getting PETSc to work, so I configured it more than just once, and in the last run I obviously forgot the option --with-scalar-type=complex.

This is what I should have done:

• Either check the log file $PETSC_DIR/arch-linux2-c-opt/conf/configure.log.

• Or take a look at the reconfigure-arch-linux2-c-opt.py.

There you can find all options you used to configure PETSc. In case you use SLEPc as well, you also need to recompile it. Now since I added the option (--with-scalar-type=complex) to the reconfigure script and ran it, everything works perfectly fine.