2

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.

Community
  • 1
  • 1
Schwartz Richard
  • 41
  • 3
  • So the dtype for `p3` is complex, right? Doesn't `PETSc` say anything about specifying a `dtype`? I'm guessing it initializes some array (`iarray`?) and copies data from `p3` to it. If it was initialized as `float`, that copy would produce the error. – hpaulj Apr 09 '15 at 02:53
  • http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscComplex.html I wonder if your PETSC is configured for COMPLEX as described on this page? – hpaulj Apr 09 '15 at 03:20
  • Thank you guys! 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. Now since I added the option to the reconfigure script and ran it, everything works perfectly fine. Sorry that I bothered you guys. – Schwartz Richard Apr 09 '15 at 09:15
  • @SchwartzRichard Could you turn that comment into an answer? – ali_m Apr 09 '15 at 13:06

1 Answers1

2

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.

ali_m
  • 71,714
  • 23
  • 223
  • 298
Schwartz Richard
  • 41
  • 3