Commits

BarryFSmith  committed e2bba1e Merge

Merge branch 'barry/add-complex-svd' into next

  • Participants
  • Parent commits d5bd065, ef47b4b
  • Branches jed/next-seaice, knepley/fix-quadrature-order 2
    1. madams/sr1
    2. next-oct-2014

Comments (0)

Files changed (1)

File src/ksp/pc/impls/svd/svd.c

     ierr = PetscFPTrapPop();CHKERRQ(ierr);
   }
 #else
-  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex");
+  {
+    PetscBLASInt lierr;
+    PetscReal    *rwork,*dd;
+    ierr = PetscMalloc(5*nb*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
+    ierr = PetscMalloc(nb*sizeof(PetscReal),&dd);CHKERRQ(ierr);
+    ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
+    PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
+    if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
+    ierr = PetscFree(rwork);CHKERRQ(ierr);
+    for (i=0; i<nb; i++) d[i] = dd[i];
+    ierr = PetscFree(dd);CHKERRQ(ierr);
+    ierr = PetscFPTrapPop();CHKERRQ(ierr);
+  }
 #endif
   ierr = MatDenseRestoreArray(jac->A,&a);CHKERRQ(ierr);
   ierr = MatDenseRestoreArray(jac->U,&u);CHKERRQ(ierr);
   PetscFunctionBegin;
   ierr = PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);CHKERRQ(ierr);
   ierr = PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);CHKERRQ(ierr);
+#if !defined(PETSC_USE_COMPLEX)
   ierr = MatMultTranspose(jac->U,xred,work);CHKERRQ(ierr);
+#else
+  ierr = MatMultHermitianTranspose(jac->U,xred,work);CHKERRQ(ierr);
+#endif
   ierr = VecPointwiseMult(work,work,jac->diag);CHKERRQ(ierr);
+#if !defined(PETSC_USE_COMPLEX)
   ierr = MatMultTranspose(jac->Vt,work,yred);CHKERRQ(ierr);
+#else
+  ierr = MatMultHermitianTranspose(jac->Vt,work,yred);CHKERRQ(ierr);
+#endif
   ierr = PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);CHKERRQ(ierr);
   ierr = PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);CHKERRQ(ierr);
   PetscFunctionReturn(0);
   PetscFunctionBegin;
   ierr = PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);CHKERRQ(ierr);
   ierr = PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);CHKERRQ(ierr);
+#if !defined(PETSC_USE_COMPLEX)
   ierr = MatMultTranspose(jac->Vt,work,yred);CHKERRQ(ierr);
+#else
+  ierr = MatMultHermitianTranspose(jac->Vt,work,yred);CHKERRQ(ierr);
+#endif
   ierr = VecPointwiseMult(work,work,jac->diag);CHKERRQ(ierr);
+#if !defined(PETSC_USE_COMPLEX)
   ierr = MatMultTranspose(jac->U,xred,work);CHKERRQ(ierr);
+#else
+  ierr = MatMultHermitianTranspose(jac->U,xred,work);CHKERRQ(ierr);
+#endif
   ierr = PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);CHKERRQ(ierr);
   ierr = PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);CHKERRQ(ierr);
   PetscFunctionReturn(0);