1 #include <../src/tao/matrix/adamat.h> /*I "mat.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "MatCreateADA" 5 /*@C 6 MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal 7 8 Collective on matrix 9 10 Input Parameters: 11 + mat - matrix of arbitrary type 12 . d1 - A vector with diagonal elements of D1 13 - d2 - A vector 14 15 Output Parameters: 16 . J - New matrix whose operations are defined in terms of mat, D1, and D2. 17 18 Notes: 19 The user provides the input data and is responsible for destroying 20 this data after matrix J has been destroyed. 21 The operation MatMult(A,D2,D1) must be well defined. 22 Before calling the operation MatGetDiagonal(), the function 23 MatADAComputeDiagonal() must be called. The matrices A and D1 must 24 be the same during calls to MatADAComputeDiagonal() and 25 MatGetDiagonal(). 26 27 Level: developer 28 29 .seealso: MatCreate() 30 @*/ 31 PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J) 32 { 33 MPI_Comm comm=((PetscObject)mat)->comm; 34 TaoMatADACtx ctx; 35 PetscErrorCode ierr; 36 PetscInt nloc,n; 37 38 PetscFunctionBegin; 39 ierr = PetscNew(&ctx);CHKERRQ(ierr); 40 ctx->A=mat; 41 ctx->D1=d1; 42 ctx->D2=d2; 43 if (d1){ 44 ierr = VecDuplicate(d1,&ctx->W);CHKERRQ(ierr); 45 ierr = PetscObjectReference((PetscObject)d1);CHKERRQ(ierr); 46 } else { 47 ctx->W = NULL; 48 } 49 if (d2){ 50 ierr = VecDuplicate(d2,&ctx->W2);CHKERRQ(ierr); 51 ierr = VecDuplicate(d2,&ctx->ADADiag);CHKERRQ(ierr); 52 ierr = PetscObjectReference((PetscObject)d2);CHKERRQ(ierr); 53 } else { 54 ctx->W2 = NULL; 55 ctx->ADADiag = NULL; 56 } 57 58 ctx->GotDiag=0; 59 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 60 61 ierr=VecGetLocalSize(d2,&nloc);CHKERRQ(ierr); 62 ierr=VecGetSize(d2,&n);CHKERRQ(ierr); 63 64 ierr = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(ierr); 65 66 ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);CHKERRQ(ierr); 67 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);CHKERRQ(ierr); 68 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);CHKERRQ(ierr); 69 ierr = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);CHKERRQ(ierr); 70 ierr = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);CHKERRQ(ierr); 71 ierr = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);CHKERRQ(ierr); 72 ierr = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);CHKERRQ(ierr); 73 ierr = MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);CHKERRQ(ierr); 74 ierr = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);CHKERRQ(ierr); 75 ierr = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);CHKERRQ(ierr); 76 ierr = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_ADA);CHKERRQ(ierr); 77 ierr = MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);CHKERRQ(ierr); 78 ierr = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);CHKERRQ(ierr); 79 ierr = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_ADA);CHKERRQ(ierr); 80 81 ierr = PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);CHKERRQ(ierr); 82 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));CHKERRQ(ierr); 83 84 ierr = MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "MatMult_ADA" 90 PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y) 91 { 92 TaoMatADACtx ctx; 93 PetscReal one = 1.0; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 98 ierr = MatMult(ctx->A,a,ctx->W);CHKERRQ(ierr); 99 if (ctx->D1){ 100 ierr = VecPointwiseMult(ctx->W,ctx->D1,ctx->W);CHKERRQ(ierr); 101 } 102 ierr = MatMultTranspose(ctx->A,ctx->W,y);CHKERRQ(ierr); 103 if (ctx->D2){ 104 ierr = VecPointwiseMult(ctx->W2, ctx->D2, a);CHKERRQ(ierr); 105 ierr = VecAXPY(y, one, ctx->W2);CHKERRQ(ierr); 106 } 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "MatMultTranspose_ADA" 112 PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y) 113 { 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 ierr = MatMult_ADA(mat,a,y);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNCT__ 122 #define __FUNCT__ "MatDiagonalSet_ADA" 123 PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M) 124 { 125 TaoMatADACtx ctx; 126 PetscReal zero=0.0,one = 1.0; 127 PetscErrorCode ierr; 128 129 PetscFunctionBegin; 130 ierr = MatShellGetContext(M,(void **)&ctx);CHKERRQ(ierr); 131 if (!ctx->D2){ 132 ierr = VecDuplicate(D,&ctx->D2);CHKERRQ(ierr); 133 ierr = VecSet(ctx->D2, zero);CHKERRQ(ierr); 134 } 135 ierr = VecAXPY(ctx->D2, one, D);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 } 138 139 #undef __FUNCT__ 140 #define __FUNCT__ "MatDestroy_ADA" 141 PetscErrorCode MatDestroy_ADA(Mat mat) 142 { 143 PetscErrorCode ierr; 144 TaoMatADACtx ctx; 145 146 PetscFunctionBegin; 147 ierr=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 148 ierr=VecDestroy(&ctx->W);CHKERRQ(ierr); 149 ierr=VecDestroy(&ctx->W2);CHKERRQ(ierr); 150 ierr=VecDestroy(&ctx->ADADiag);CHKERRQ(ierr); 151 ierr=MatDestroy(&ctx->A);CHKERRQ(ierr); 152 ierr=VecDestroy(&ctx->D1);CHKERRQ(ierr); 153 ierr=VecDestroy(&ctx->D2);CHKERRQ(ierr); 154 ierr = PetscFree(ctx);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "MatView_ADA" 160 PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer) 161 { 162 PetscFunctionBegin; 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "MatShift_ADA" 168 PetscErrorCode MatShift_ADA(Mat Y, PetscReal a) 169 { 170 PetscErrorCode ierr; 171 TaoMatADACtx ctx; 172 173 PetscFunctionBegin; 174 ierr = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(ierr); 175 ierr = VecShift(ctx->D2,a);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 179 #undef __FUNCT__ 180 #define __FUNCT__ "MatDuplicate_ADA" 181 PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M) 182 { 183 PetscErrorCode ierr; 184 TaoMatADACtx ctx; 185 Mat A2; 186 Vec D1b=NULL,D2b; 187 188 PetscFunctionBegin; 189 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 190 ierr = MatDuplicate(ctx->A,op,&A2);CHKERRQ(ierr); 191 if (ctx->D1){ 192 ierr = VecDuplicate(ctx->D1,&D1b);CHKERRQ(ierr); 193 ierr = VecCopy(ctx->D1,D1b);CHKERRQ(ierr); 194 } 195 ierr = VecDuplicate(ctx->D2,&D2b);CHKERRQ(ierr); 196 ierr = VecCopy(ctx->D2,D2b);CHKERRQ(ierr); 197 ierr = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(ierr); 198 if (ctx->D1){ 199 ierr = PetscObjectDereference((PetscObject)D1b);CHKERRQ(ierr); 200 } 201 ierr = PetscObjectDereference((PetscObject)D2b);CHKERRQ(ierr); 202 ierr = PetscObjectDereference((PetscObject)A2);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "MatEqual_ADA" 208 PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg) 209 { 210 PetscErrorCode ierr; 211 TaoMatADACtx ctx1,ctx2; 212 213 PetscFunctionBegin; 214 ierr = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(ierr); 215 ierr = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(ierr); 216 ierr = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(ierr); 217 if (*flg==PETSC_TRUE){ 218 ierr = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(ierr); 219 } 220 if (*flg==PETSC_TRUE){ 221 ierr = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(ierr); 222 } 223 PetscFunctionReturn(0); 224 } 225 226 #undef __FUNCT__ 227 #define __FUNCT__ "MatScale_ADA" 228 PetscErrorCode MatScale_ADA(Mat mat, PetscReal a) 229 { 230 PetscErrorCode ierr; 231 TaoMatADACtx ctx; 232 233 PetscFunctionBegin; 234 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 235 ierr = VecScale(ctx->D1,a);CHKERRQ(ierr); 236 if (ctx->D2){ 237 ierr = VecScale(ctx->D2,a);CHKERRQ(ierr); 238 } 239 PetscFunctionReturn(0); 240 } 241 242 #undef __FUNCT__ 243 #define __FUNCT__ "MatTranspose_ADA" 244 PetscErrorCode MatTranspose_ADA(Mat mat,Mat *B) 245 { 246 PetscErrorCode ierr; 247 TaoMatADACtx ctx; 248 249 PetscFunctionBegin; 250 if (*B){ 251 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 252 ierr = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(ierr); 253 } 254 PetscFunctionReturn(0); 255 } 256 257 #undef __FUNCT__ 258 #define __FUNCT__ "MatADAComputeDiagonal" 259 PetscErrorCode MatADAComputeDiagonal(Mat mat) 260 { 261 PetscErrorCode ierr; 262 PetscInt i,m,n,low,high; 263 PetscScalar *dtemp,*dptr; 264 TaoMatADACtx ctx; 265 266 PetscFunctionBegin; 267 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 268 ierr = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(ierr); 269 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 270 271 ierr = PetscMalloc1(n,&dtemp );CHKERRQ(ierr); 272 273 for (i=0; i<n; i++){ 274 ierr = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(ierr); 275 ierr = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(ierr); 276 ierr = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr); 277 } 278 for (i=0; i<n; i++){ 279 ierr = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr); 280 } 281 282 ierr = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(ierr); 283 for (i=low; i<high; i++){ 284 dptr[i-low]= dtemp[i]; 285 } 286 ierr = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(ierr); 287 ierr = PetscFree(dtemp);CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatGetDiagonal_ADA" 293 PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v) 294 { 295 PetscErrorCode ierr; 296 PetscReal one=1.0; 297 TaoMatADACtx ctx; 298 299 PetscFunctionBegin; 300 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 301 ierr = MatADAComputeDiagonal(mat);CHKERRQ(ierr); 302 ierr=VecCopy(ctx->ADADiag,v);CHKERRQ(ierr); 303 if (ctx->D2){ 304 ierr=VecAXPY(v, one, ctx->D2);CHKERRQ(ierr); 305 } 306 PetscFunctionReturn(0); 307 } 308 309 #undef __FUNCT__ 310 #define __FUNCT__ "MatGetSubMatrices_ADA" 311 PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B) 312 { 313 PetscErrorCode ierr; 314 PetscInt i; 315 316 PetscFunctionBegin; 317 if (scall == MAT_INITIAL_MATRIX) { 318 ierr = PetscMalloc1(n+1,B );CHKERRQ(ierr); 319 } 320 for ( i=0; i<n; i++ ) { 321 ierr = MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 322 } 323 PetscFunctionReturn(0); 324 } 325 326 #undef __FUNCT__ 327 #define __FUNCT__ "MatGetSubMatrix_ADA" 328 PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat) 329 { 330 PetscErrorCode ierr; 331 PetscInt low,high; 332 PetscInt n,nlocal,i; 333 const PetscInt *iptr; 334 const PetscScalar *dptr; 335 PetscScalar *ddptr,zero=0.0; 336 VecType type_name; 337 IS ISrow; 338 Vec D1,D2; 339 Mat Atemp; 340 TaoMatADACtx ctx; 341 342 PetscFunctionBegin; 343 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 344 345 ierr = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(ierr); 346 ierr = ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);CHKERRQ(ierr); 347 ierr = MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);CHKERRQ(ierr); 348 ierr = ISDestroy(&ISrow);CHKERRQ(ierr); 349 350 if (ctx->D1){ 351 ierr=VecDuplicate(ctx->D1,&D1);CHKERRQ(ierr); 352 ierr=VecCopy(ctx->D1,D1);CHKERRQ(ierr); 353 } else { 354 D1 = NULL; 355 } 356 357 if (ctx->D2){ 358 ierr=ISGetSize(isrow,&n);CHKERRQ(ierr); 359 ierr=ISGetLocalSize(isrow,&nlocal);CHKERRQ(ierr); 360 ierr=VecCreate(((PetscObject)(ctx->D2))->comm,&D2);CHKERRQ(ierr); 361 ierr=VecGetType(ctx->D2,&type_name);CHKERRQ(ierr); 362 ierr=VecSetSizes(D2,nlocal,n);CHKERRQ(ierr); 363 ierr=VecSetType(D2,type_name);CHKERRQ(ierr); 364 ierr=VecSet(D2, zero);CHKERRQ(ierr); 365 ierr=VecGetArrayRead(ctx->D2, &dptr);CHKERRQ(ierr); 366 ierr=VecGetArray(D2, &ddptr);CHKERRQ(ierr); 367 ierr=ISGetIndices(isrow,&iptr);CHKERRQ(ierr); 368 for (i=0;i<nlocal;i++){ 369 ddptr[i] = dptr[iptr[i]-low]; 370 } 371 ierr=ISRestoreIndices(isrow,&iptr);CHKERRQ(ierr); 372 ierr=VecRestoreArray(D2, &ddptr);CHKERRQ(ierr); 373 ierr=VecRestoreArrayRead(ctx->D2, &dptr);CHKERRQ(ierr); 374 } else { 375 D2 = NULL; 376 } 377 378 ierr = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(ierr); 379 ierr = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(ierr); 380 ierr = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(ierr); 381 if (ctx->D1){ 382 ierr = PetscObjectDereference((PetscObject)D1);CHKERRQ(ierr); 383 } 384 if (ctx->D2){ 385 ierr = PetscObjectDereference((PetscObject)D2);CHKERRQ(ierr); 386 } 387 PetscFunctionReturn(0); 388 } 389 390 #undef __FUNCT__ 391 #define __FUNCT__ "MatGetRowADA" 392 PetscErrorCode MatGetRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals) 393 { 394 PetscErrorCode ierr; 395 PetscInt m,n; 396 397 PetscFunctionBegin; 398 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 399 400 if (*ncols>0){ 401 ierr = PetscMalloc1(*ncols,cols );CHKERRQ(ierr); 402 ierr = PetscMalloc1(*ncols,vals );CHKERRQ(ierr); 403 } else { 404 *cols=NULL; 405 *vals=NULL; 406 } 407 PetscFunctionReturn(0); 408 } 409 410 #undef __FUNCT__ 411 #define __FUNCT__ "MatRestoreRowADA" 412 PetscErrorCode MatRestoreRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals) 413 { 414 PetscErrorCode ierr; 415 416 PetscFunctionBegin; 417 if (*ncols>0){ 418 ierr = PetscFree(*cols); CHKERRQ(ierr); 419 ierr = PetscFree(*vals); CHKERRQ(ierr); 420 } 421 *cols=NULL; 422 *vals=NULL; 423 PetscFunctionReturn(0); 424 } 425 426 #undef __FUNCT__ 427 #define __FUNCT__ "MatGetColumnVectorADA" 428 PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col) 429 { 430 PetscErrorCode ierr; 431 PetscInt low,high; 432 PetscScalar zero=0.0,one=1.0; 433 434 PetscFunctionBegin; 435 ierr=VecSet(Y, zero);CHKERRQ(ierr); 436 ierr=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(ierr); 437 if (col>=low && col<high){ 438 ierr=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(ierr); 439 } 440 ierr=VecAssemblyBegin(Y);CHKERRQ(ierr); 441 ierr=VecAssemblyEnd(Y);CHKERRQ(ierr); 442 ierr=MatMult_ADA(mat,Y,Y);CHKERRQ(ierr); 443 PetscFunctionReturn(0); 444 } 445 446 PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat) 447 { 448 PetscErrorCode ierr; 449 PetscMPIInt size; 450 PetscBool sametype, issame, isdense, isseqdense; 451 TaoMatADACtx ctx; 452 453 PetscFunctionBegin; 454 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 455 ierr = MPI_Comm_size(((PetscObject)mat)->comm,&size);CHKERRQ(ierr); 456 457 ierr = PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr); 458 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);CHKERRQ(ierr); 459 ierr = PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);CHKERRQ(ierr); 460 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 461 462 if (sametype || issame) { 463 ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(ierr); 464 } else if (isdense) { 465 PetscInt i,j,low,high,m,n,M,N; 466 const PetscScalar *dptr; 467 Vec X; 468 469 ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr); 470 ierr=MatGetSize(mat,&M,&N);CHKERRQ(ierr); 471 ierr=MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 472 ierr = MatCreateDense(((PetscObject)mat)->comm,m,m,N,N,NULL,NewMat);CHKERRQ(ierr); 473 ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr); 474 for (i=0;i<M;i++){ 475 ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr); 476 ierr = VecGetArrayRead(X,&dptr);CHKERRQ(ierr); 477 for (j=0; j<high-low; j++){ 478 ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr); 479 } 480 ierr = VecRestoreArrayRead(X,&dptr);CHKERRQ(ierr); 481 } 482 ierr = MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 483 ierr = MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 484 ierr = VecDestroy(&X);CHKERRQ(ierr); 485 } else if (isseqdense && size==1){ 486 PetscInt i,j,low,high,m,n,M,N; 487 const PetscScalar *dptr; 488 Vec X; 489 490 ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr); 491 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 492 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 493 ierr = MatCreateSeqDense(((PetscObject)mat)->comm,N,N,NULL,NewMat);CHKERRQ(ierr); 494 ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr); 495 for (i=0;i<M;i++){ 496 ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr); 497 ierr = VecGetArrayRead(X,&dptr);CHKERRQ(ierr); 498 for (j=0; j<high-low; j++){ 499 ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr); 500 } 501 ierr = VecRestoreArrayRead(X,&dptr);CHKERRQ(ierr); 502 } 503 ierr = MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 504 ierr = MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 505 ierr = VecDestroy(&X);CHKERRQ(ierr); 506 } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type"); 507 PetscFunctionReturn(0); 508 } 509 510 #undef __FUNCT__ 511 #define __FUNCT__ "MatNorm_ADA" 512 PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm) 513 { 514 PetscErrorCode ierr; 515 TaoMatADACtx ctx; 516 517 PetscFunctionBegin; 518 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 519 if (type == NORM_FROBENIUS) { 520 *norm = 1.0; 521 } else if (type == NORM_1 || type == NORM_INFINITY) { 522 *norm = 1.0; 523 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 524 PetscFunctionReturn(0); 525 } 526