xref: /petsc/src/mat/impls/lrc/lrc.c (revision ceb9aaf78cdd026d60598f6b38b8f23191eb747b) 
1*ceb9aaf7SBarry Smith #define PETSCMAT_DLL
2*ceb9aaf7SBarry Smith 
3*ceb9aaf7SBarry Smith #include "src/mat/matimpl.h"          /*I "petscmat.h" I*/
4*ceb9aaf7SBarry Smith 
5*ceb9aaf7SBarry Smith typedef struct {
6*ceb9aaf7SBarry Smith   Mat A;
7*ceb9aaf7SBarry Smith   Vec w;
8*ceb9aaf7SBarry Smith } Mat_Normal;
9*ceb9aaf7SBarry Smith 
10*ceb9aaf7SBarry Smith #undef __FUNCT__
11*ceb9aaf7SBarry Smith #define __FUNCT__ "MatMult_Normal"
12*ceb9aaf7SBarry Smith PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y)
13*ceb9aaf7SBarry Smith {
14*ceb9aaf7SBarry Smith   Mat_Normal     *Na = (Mat_Normal*)N->data;
15*ceb9aaf7SBarry Smith   PetscErrorCode ierr;
16*ceb9aaf7SBarry Smith 
17*ceb9aaf7SBarry Smith   PetscFunctionBegin;
18*ceb9aaf7SBarry Smith   ierr = MatMult(Na->A,x,Na->w);CHKERRQ(ierr);
19*ceb9aaf7SBarry Smith   ierr = MatMultTranspose(Na->A,Na->w,y);CHKERRQ(ierr);
20*ceb9aaf7SBarry Smith   PetscFunctionReturn(0);
21*ceb9aaf7SBarry Smith }
22*ceb9aaf7SBarry Smith 
23*ceb9aaf7SBarry Smith #undef __FUNCT__
24*ceb9aaf7SBarry Smith #define __FUNCT__ "MatDestroy_Normal"
25*ceb9aaf7SBarry Smith PetscErrorCode MatDestroy_Normal(Mat N)
26*ceb9aaf7SBarry Smith {
27*ceb9aaf7SBarry Smith   Mat_Normal     *Na = (Mat_Normal*)N->data;
28*ceb9aaf7SBarry Smith   PetscErrorCode ierr;
29*ceb9aaf7SBarry Smith 
30*ceb9aaf7SBarry Smith   PetscFunctionBegin;
31*ceb9aaf7SBarry Smith   ierr = PetscObjectDereference((PetscObject)Na->A);CHKERRQ(ierr);
32*ceb9aaf7SBarry Smith   ierr = VecDestroy(Na->w);CHKERRQ(ierr);
33*ceb9aaf7SBarry Smith   ierr = PetscFree(Na);CHKERRQ(ierr);
34*ceb9aaf7SBarry Smith   PetscFunctionReturn(0);
35*ceb9aaf7SBarry Smith }
36*ceb9aaf7SBarry Smith 
37*ceb9aaf7SBarry Smith /*
38*ceb9aaf7SBarry Smith       Slow, nonscalable version
39*ceb9aaf7SBarry Smith */
40*ceb9aaf7SBarry Smith #undef __FUNCT__
41*ceb9aaf7SBarry Smith #define __FUNCT__ "MatGetDiagonal_Normal"
42*ceb9aaf7SBarry Smith PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v)
43*ceb9aaf7SBarry Smith {
44*ceb9aaf7SBarry Smith   Mat_Normal        *Na = (Mat_Normal*)N->data;
45*ceb9aaf7SBarry Smith   Mat               A = Na->A;
46*ceb9aaf7SBarry Smith   PetscErrorCode    ierr;
47*ceb9aaf7SBarry Smith   PetscInt          i,j,rstart,rend,nnz;
48*ceb9aaf7SBarry Smith   const PetscInt    *cols;
49*ceb9aaf7SBarry Smith   PetscScalar       *diag,*work,*values;
50*ceb9aaf7SBarry Smith   const PetscScalar *mvalues;
51*ceb9aaf7SBarry Smith   PetscMap          cmap;
52*ceb9aaf7SBarry Smith 
53*ceb9aaf7SBarry Smith   PetscFunctionBegin;
54*ceb9aaf7SBarry Smith   ierr = PetscMalloc(2*A->N*sizeof(PetscScalar),&diag);CHKERRQ(ierr);
55*ceb9aaf7SBarry Smith   work = diag + A->N;
56*ceb9aaf7SBarry Smith   ierr = PetscMemzero(work,A->N*sizeof(PetscScalar));CHKERRQ(ierr);
57*ceb9aaf7SBarry Smith   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
58*ceb9aaf7SBarry Smith   for (i=rstart; i<rend; i++) {
59*ceb9aaf7SBarry Smith     ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
60*ceb9aaf7SBarry Smith     for (j=0; j<nnz; j++) {
61*ceb9aaf7SBarry Smith       work[cols[j]] += mvalues[j]*mvalues[j];
62*ceb9aaf7SBarry Smith     }
63*ceb9aaf7SBarry Smith     ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
64*ceb9aaf7SBarry Smith   }
65*ceb9aaf7SBarry Smith   ierr = MPI_Allreduce(work,diag,A->N,MPIU_SCALAR,MPI_SUM,N->comm);CHKERRQ(ierr);
66*ceb9aaf7SBarry Smith   ierr = MatGetPetscMaps(A,PETSC_NULL,&cmap);CHKERRQ(ierr);
67*ceb9aaf7SBarry Smith   ierr = PetscMapGetLocalRange(cmap,&rstart,&rend);CHKERRQ(ierr);
68*ceb9aaf7SBarry Smith   ierr = VecGetArray(v,&values);CHKERRQ(ierr);
69*ceb9aaf7SBarry Smith   ierr = PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));CHKERRQ(ierr);
70*ceb9aaf7SBarry Smith   ierr = VecRestoreArray(v,&values);CHKERRQ(ierr);
71*ceb9aaf7SBarry Smith   ierr = PetscFree(diag);CHKERRQ(ierr);
72*ceb9aaf7SBarry Smith   PetscFunctionReturn(0);
73*ceb9aaf7SBarry Smith }
74*ceb9aaf7SBarry Smith 
75*ceb9aaf7SBarry Smith #undef __FUNCT__
76*ceb9aaf7SBarry Smith #define __FUNCT__ "MatCreateNormal"
77*ceb9aaf7SBarry Smith /*@
78*ceb9aaf7SBarry Smith       MatCreateNormal - Creates a new matrix object that behaves like A'*A.
79*ceb9aaf7SBarry Smith 
80*ceb9aaf7SBarry Smith    Collective on Mat
81*ceb9aaf7SBarry Smith 
82*ceb9aaf7SBarry Smith    Input Parameter:
83*ceb9aaf7SBarry Smith .   A  - the (possibly rectangular) matrix
84*ceb9aaf7SBarry Smith 
85*ceb9aaf7SBarry Smith    Output Parameter:
86*ceb9aaf7SBarry Smith .   N - the matrix that represents A'*A
87*ceb9aaf7SBarry Smith 
88*ceb9aaf7SBarry Smith    Level: intermediate
89*ceb9aaf7SBarry Smith 
90*ceb9aaf7SBarry Smith    Notes: The product A'*A is NOT actually formed! Rather the new matrix
91*ceb9aaf7SBarry Smith           object performs the matrix-vector product by first multiplying by
92*ceb9aaf7SBarry Smith           A and then A'
93*ceb9aaf7SBarry Smith @*/
94*ceb9aaf7SBarry Smith PetscErrorCode PETSCMAT_DLLEXPORT MatCreateNormal(Mat A,Mat *N)
95*ceb9aaf7SBarry Smith {
96*ceb9aaf7SBarry Smith   PetscErrorCode ierr;
97*ceb9aaf7SBarry Smith   PetscInt       m,n;
98*ceb9aaf7SBarry Smith   Mat_Normal     *Na;
99*ceb9aaf7SBarry Smith 
100*ceb9aaf7SBarry Smith   PetscFunctionBegin;
101*ceb9aaf7SBarry Smith   ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
102*ceb9aaf7SBarry Smith   ierr = MatCreate(A->comm,N);CHKERRQ(ierr);
103*ceb9aaf7SBarry Smith   ierr = MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
104*ceb9aaf7SBarry Smith   ierr = PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);CHKERRQ(ierr);
105*ceb9aaf7SBarry Smith 
106*ceb9aaf7SBarry Smith   ierr      = PetscNew(Mat_Normal,&Na);CHKERRQ(ierr);
107*ceb9aaf7SBarry Smith   Na->A     = A;
108*ceb9aaf7SBarry Smith   ierr      = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
109*ceb9aaf7SBarry Smith   (*N)->data = (void*) Na;
110*ceb9aaf7SBarry Smith 
111*ceb9aaf7SBarry Smith   ierr    = VecCreateMPI(A->comm,m,PETSC_DECIDE,&Na->w);CHKERRQ(ierr);
112*ceb9aaf7SBarry Smith   (*N)->ops->destroy     = MatDestroy_Normal;
113*ceb9aaf7SBarry Smith   (*N)->ops->mult        = MatMult_Normal;
114*ceb9aaf7SBarry Smith   (*N)->ops->getdiagonal = MatGetDiagonal_Normal;
115*ceb9aaf7SBarry Smith   (*N)->assembled        = PETSC_TRUE;
116*ceb9aaf7SBarry Smith   (*N)->N                = A->N;
117*ceb9aaf7SBarry Smith   (*N)->M                = A->N;
118*ceb9aaf7SBarry Smith   (*N)->n                = A->n;
119*ceb9aaf7SBarry Smith   (*N)->m                = A->n;
120*ceb9aaf7SBarry Smith   PetscFunctionReturn(0);
121*ceb9aaf7SBarry Smith }
122*ceb9aaf7SBarry Smith 
123