normm.c (d8e47b638cf8f604a99e9678e1df24f82d959cd7) normm.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_Normal;
8

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

139 Slow, nonscalable version
140*/
141static PetscErrorCode MatGetDiagonal_Normal(Mat N, Vec v)
142{
143 Mat_Normal *Na;
144 Mat A;
145 PetscInt i, j, rstart, rend, nnz;
146 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_Normal;
8

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

139 Slow, nonscalable version
140*/
141static PetscErrorCode MatGetDiagonal_Normal(Mat N, Vec v)
142{
143 Mat_Normal *Na;
144 Mat A;
145 PetscInt i, j, rstart, rend, nnz;
146 const PetscInt *cols;
147 PetscScalar *diag, *work, *values;
147 PetscScalar *work, *values;
148 const PetscScalar *mvalues;
148 const PetscScalar *mvalues;
149 PetscMPIInt iN;
150
151 PetscFunctionBegin;
152 PetscCall(MatShellGetContext(N, &Na));
153 A = Na->A;
149
150 PetscFunctionBegin;
151 PetscCall(MatShellGetContext(N, &Na));
152 A = Na->A;
154 PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
153 PetscCall(PetscMalloc1(A->cmap->N, &work));
155 PetscCall(PetscArrayzero(work, A->cmap->N));
156 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
157 for (i = rstart; i < rend; i++) {
158 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
159 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * mvalues[j];
160 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
161 }
154 PetscCall(PetscArrayzero(work, A->cmap->N));
155 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
156 for (i = rstart; i < rend; i++) {
157 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
158 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * mvalues[j];
159 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
160 }
162 PetscCall(PetscMPIIntCast(A->cmap->N, &iN));
163 PetscCallMPI(MPIU_Allreduce(work, diag, iN, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
161 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, work, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
164 rstart = N->cmap->rstart;
165 rend = N->cmap->rend;
166 PetscCall(VecGetArray(v, &values));
162 rstart = N->cmap->rstart;
163 rend = N->cmap->rend;
164 PetscCall(VecGetArray(v, &values));
167 PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
165 PetscCall(PetscArraycpy(values, work + rstart, rend - rstart));
168 PetscCall(VecRestoreArray(v, &values));
166 PetscCall(VecRestoreArray(v, &values));
169 PetscCall(PetscFree2(diag, work));
167 PetscCall(PetscFree(work));
170 PetscFunctionReturn(PETSC_SUCCESS);
171}
172
173static PetscErrorCode MatGetDiagonalBlock_Normal(Mat N, Mat *D)
174{
175 Mat_Normal *Na;
176 Mat M, A;
177

--- 270 unchanged lines hidden --- 		168 PetscFunctionReturn(PETSC_SUCCESS);
169}
170
171static PetscErrorCode MatGetDiagonalBlock_Normal(Mat N, Mat *D)
172{
173 Mat_Normal *Na;
174 Mat M, A;
175

--- 270 unchanged lines hidden ---