normmh.c (d8e47b638cf8f604a99e9678e1df24f82d959cd7) normmh.c (e91c04dfc8a52dee1965211bb1cc8e5bf775178f)
1#include <../src/mat/impls/shell/shell.h> /*I "petscmat.h" I*/
2
3typedef struct {
4 Mat A;
5 Mat D; /* local submatrix for diagonal part */
6 Vec w;
7} Mat_NormalHermitian;
8

--- 114 unchanged lines hidden (view full) ---

123 Slow, nonscalable version
124*/
125static PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
126{
127 Mat_NormalHermitian *Na;
128 Mat A;
129 PetscInt i, j, rstart, rend, nnz;
130 const PetscInt *cols;
1#include <../src/mat/impls/shell/shell.h> /*I "petscmat.h" I*/
2
3typedef struct {
4 Mat A;
5 Mat D; /* local submatrix for diagonal part */
6 Vec w;
7} Mat_NormalHermitian;
8

--- 114 unchanged lines hidden (view full) ---

123 Slow, nonscalable version
124*/
125static PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
126{
127 Mat_NormalHermitian *Na;
128 Mat A;
129 PetscInt i, j, rstart, rend, nnz;
130 const PetscInt *cols;
131 PetscScalar *diag, *work, *values;
131 PetscScalar *work, *values;
132 const PetscScalar *mvalues;
132 const PetscScalar *mvalues;
133 PetscMPIInt iN;
134
135 PetscFunctionBegin;
136 PetscCall(MatShellGetContext(N, &Na));
137 A = Na->A;
133
134 PetscFunctionBegin;
135 PetscCall(MatShellGetContext(N, &Na));
136 A = Na->A;
138 PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
137 PetscCall(PetscMalloc1(A->cmap->N, &work));
139 PetscCall(PetscArrayzero(work, A->cmap->N));
140 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
141 for (i = rstart; i < rend; i++) {
142 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
143 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
144 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
145 }
138 PetscCall(PetscArrayzero(work, A->cmap->N));
139 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
140 for (i = rstart; i < rend; i++) {
141 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
142 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
143 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
144 }
146 PetscCall(PetscMPIIntCast(A->cmap->N, &iN));
147 PetscCallMPI(MPIU_Allreduce(work, diag, iN, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
145 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, work, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
148 rstart = N->cmap->rstart;
149 rend = N->cmap->rend;
150 PetscCall(VecGetArray(v, &values));
146 rstart = N->cmap->rstart;
147 rend = N->cmap->rend;
148 PetscCall(VecGetArray(v, &values));
151 PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
149 PetscCall(PetscArraycpy(values, work + rstart, rend - rstart));
152 PetscCall(VecRestoreArray(v, &values));
150 PetscCall(VecRestoreArray(v, &values));
153 PetscCall(PetscFree2(diag, work));
151 PetscCall(PetscFree(work));
154 PetscFunctionReturn(PETSC_SUCCESS);
155}
156
157static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
158{
159 Mat_NormalHermitian *Na;
160 Mat M, A;
161

--- 175 unchanged lines hidden --- 		152 PetscFunctionReturn(PETSC_SUCCESS);
153}
154
155static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
156{
157 Mat_NormalHermitian *Na;
158 Mat M, A;
159

--- 175 unchanged lines hidden ---