<html><head><meta http-equiv="content-type" content="text/html; charset=utf-8"></head><body style="overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div><br></div> It is simplier than you are making it out to be. Each IS[] is a list of rows (and columns) in the sub (domain) matrix. In your case with the matrix of 144 by 144 the indices will go from 0 to 143.<div><br></div><div> In your simple Fortran code you have a completely different problem. A matrix with 8 rows and columns. In that case if you want the first IS to represent just the first row (and column) in the matrix then it should contain only 0. The second submatrix which is all rows (but the first) should have 1,2,3,4,5,6,7</div><div><br></div><div> I do not understand why your code has </div><div><br></div><div><blockquote type="cite"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div dir="ltr"><div class="gmail_quote"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div dir="ltr"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div dir="ltr"><font face="monospace">indices_first_domain = [0,1,8,9] ! corresponds to diag(1,1)</font></div></blockquote></div></div></blockquote></div></div></div></blockquote></div></blockquote><br></div><div>it should just be 0</div><div><br></div><div><br><div><br></div><div><br><div><br><blockquote type="cite"><div>On May 9, 2023, at 12:44 PM, LEONARDO MUTTI <leonardo.mutti01@universitadipavia.it> wrote:</div><br class="Apple-interchange-newline"><div><div dir="ltr">Partial typo: I expect 9x(16+16) numbers to be stored in subdomain_IS : # subdomains x (row indices of the submatrix + col indices of the submatrix).</div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno mar 9 mag 2023 alle ore 18:31 LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it">leonardo.mutti01@universitadipavia.it</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">---------- Forwarded message ---------<br>Da: <strong class="gmail_sendername" dir="auto">LEONARDO MUTTI</strong> <span dir="auto"><<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>></span><br>Date: mar 9 mag 2023 alle ore 18:29<br>Subject: Re: [petsc-users] Understanding index sets for PCGASM<br>To: Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br></div><br><br><div dir="ltr">Thank you for your answer, but I am still confused, sorry. <div>Consider <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90" target="_blank">https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90</a> on one processor.</div><div>Let M=12 for the sake of simplicity, i.e. we deal with a 12x12 2D grid, hence, a 144x144 matrix.</div><div>Let NSubx = 3, so that on the grid we do 3 vertical and 3 horizontal subdivisions.</div><div>We should obtain 9 subdomains that are grids of 4x4 nodes each, thus corresponding to 9 submatrices of size 16x16. </div><div><div>In my run I obtain NSub = 9 (great) and subdomain_IS(i), i=1,...,9, reads:</div></div><div><br></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><i>IS Object: 1 MPI process</i><br></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 0</i></div><div><i>1 1</i></div><div><i>2 2</i></div><div><i>3 3</i></div><div><i>4 12</i></div><div><i>5 13</i></div><div><i>6 14</i></div><div><i>7 15</i></div><div><i>8 24</i></div><div><i>9 25</i></div><div><i>10 26</i></div><div><i>11 27</i></div><div><i>12 36</i></div><div><i>13 37</i></div><div><i>14 38</i></div><div><i>15 39</i></div><div><i>IS Object: 1 MPI process</i></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 4</i></div><div><i>1 5</i></div><div><i>2 6</i></div><div><i>3 7</i></div><div><i>4 16</i></div><div><i>5 17</i></div><div><i>6 18</i></div><div><i>7 19</i></div><div><i>8 28</i></div><div><i>9 29</i></div><div><i>10 30</i></div><div><i>11 31</i></div><div><i>12 40</i></div><div><i>13 41</i></div><div><i>14 42</i></div><div><i>15 43</i></div><div><i>IS Object: 1 MPI process</i></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 8</i></div><div><i>1 9</i></div><div><i>2 10</i></div><div><i>3 11</i></div><div><i>4 20</i></div><div><i>5 21</i></div><div><i>6 22</i></div><div><i>7 23</i></div><div><i>8 32</i></div><div><i>9 33</i></div><div><i>10 34</i></div><div><i>11 35</i></div><div><i>12 44</i></div><div><i>13 45</i></div><div><i>14 46</i></div><div><i>15 47</i></div><div><i>IS Object: 1 MPI process</i></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 48</i></div><div><i>1 49</i></div><div><i>2 50</i></div><div><i>3 51</i></div><div><i>4 60</i></div><div><i>5 61</i></div><div><i>6 62</i></div><div><i>7 63</i></div><div><i>8 72</i></div><div><i>9 73</i></div><div><i>10 74</i></div><div><i>11 75</i></div><div><i>12 84</i></div><div><i>13 85</i></div><div><i>14 86</i></div><div><i>15 87</i></div><div><i>IS Object: 1 MPI process</i></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 52</i></div><div><i>1 53</i></div><div><i>2 54</i></div><div><i>3 55</i></div><div><i>4 64</i></div><div><i>5 65</i></div><div><i>6 66</i></div><div><i>7 67</i></div><div><i>8 76</i></div><div><i>9 77</i></div><div><i>10 78</i></div><div><i>11 79</i></div><div><i>12 88</i></div><div><i>13 89</i></div><div><i>14 90</i></div><div><i>15 91</i></div><div><i>IS Object: 1 MPI process</i></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 56</i></div><div><i>1 57</i></div><div><i>2 58</i></div><div><i>3 59</i></div><div><i>4 68</i></div><div><i>5 69</i></div><div><i>6 70</i></div><div><i>7 71</i></div><div><i>8 80</i></div><div><i>9 81</i></div><div><i>10 82</i></div><div><i>11 83</i></div><div><i>12 92</i></div><div><i>13 93</i></div><div><i>14 94</i></div><div><i>15 95</i></div><div><i>IS Object: 1 MPI process</i></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 96</i></div><div><i>1 97</i></div><div><i>2 98</i></div><div><i>3 99</i></div><div><i>4 108</i></div><div><i>5 109</i></div><div><i>6 110</i></div><div><i>7 111</i></div><div><i>8 120</i></div><div><i>9 121</i></div><div><i>10 122</i></div><div><i>11 123</i></div><div><i>12 132</i></div><div><i>13 133</i></div><div><i>14 134</i></div><div><i>15 135</i></div><div><i>IS Object: 1 MPI process</i></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 100</i></div><div><i>1 101</i></div><div><i>2 102</i></div><div><i>3 103</i></div><div><i>4 112</i></div><div><i>5 113</i></div><div><i>6 114</i></div><div><i>7 115</i></div><div><i>8 124</i></div><div><i>9 125</i></div><div><i>10 126</i></div><div><i>11 127</i></div><div><i>12 136</i></div><div><i>13 137</i></div><div><i>14 138</i></div><div><i>15 139</i></div><div><i>IS Object: 1 MPI process</i></div><div><i> type: general</i></div><div><i>Number of indices in set 16</i></div><div><i>0 104</i></div><div><i>1 105</i></div><div><i>2 106</i></div><div><i>3 107</i></div><div><i>4 116</i></div><div><i>5 117</i></div><div><i>6 118</i></div><div><i>7 119</i></div><div><i>8 128</i></div><div><i>9 129</i></div><div><i>10 130</i></div><div><i>11 131</i></div><div><i>12 140</i></div><div><i>13 141</i></div><div><i>14 142</i></div><div><i>15 143</i></div><div><br></div></blockquote><div>As you said, no number here reaches 144. </div><div>But the number stored in subdomain_IS are 9x16= #subdomains x 16, whereas I would expect, also given your latest reply, 9x16x16x2=#subdomains x submatrix height x submatrix width x length of a (row,column) pair.</div><div>It would really help me if you could briefly explain how the output above encodes the subdivision into subdomains.<br></div><div>Many thanks again,<br></div><div>Leonardo<br><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><br></div></blockquote></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno mar 9 mag 2023 alle ore 16:24 Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Tue, May 9, 2023 at 10:05 AM LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Great thanks! I can now successfully run <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90" target="_blank">https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90</a>.</div><div><br>Going forward with my experiments, let me post a new code snippet (very similar to
ex71f.F90) that I cannot get to work, probably I must be setting up the IS objects incorrectly.</div><div><br></div><div>I have an 8x8 matrix A=diag(1,1,2,2,...,2) and a vector b=(0.5,...,0.5). We have only one processor, and I want to solve Ax=b using GASM. In particular, KSP is set to preonly, GASM is the preconditioner and it uses on each submatrix an lu direct solver (sub_ksp = preonly, sub_pc = lu).</div><div><br></div><div>For the GASM algorithm, I divide A into diag(1,1) and diag(2,2,...,2). For simplicity I set 0 overlap. Now I want to use GASM to solve Ax=b. The code follows.</div><div><font face="monospace"><br>#include <petsc/finclude/petscmat.h><br>#include <petsc/finclude/petscksp.h><br>#include <petsc/finclude/petscpc.h><br> USE petscmat<br> USE petscksp<br> USE petscpc<br> USE MPI<br><br></font></div><div><font face="monospace"> Mat :: A<br> Vec :: b, x<br> PetscInt :: M, I, J, ISLen, NSub<br> PetscMPIInt :: size<br> PetscErrorCode :: ierr<br> PetscScalar :: v<br> KSP :: ksp<br> PC :: pc<br> IS :: subdomains_IS(2), inflated_IS(2)<br> PetscInt,DIMENSION(4) :: indices_first_domain<br> PetscInt,DIMENSION(36) :: indices_second_domain<br> <br> call PetscInitialize(PETSC_NULL_CHARACTER, ierr)<br> call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)<br><br><br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! INTRO: create matrix and right hand side<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> <br> WRITE(*,*) "Assembling A,b"<br><br> M = 8<br> call MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,<br> & M, M, PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,<br> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,A, ierr)<br> DO I=1,M<br> DO J=1,M<br> IF ((I .EQ. J) .AND. (I .LE. 2 )) THEN<br> v = 1<br> ELSE IF ((I .EQ. J) .AND. (I .GT. 2 )) THEN<br> v = 2<br> ELSE<br> v = 0<br> ENDIF<br> call MatSetValue(A, I-1, J-1, v, INSERT_VALUES, ierr)<br> END DO<br> END DO<br><br> call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)<br> call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)<br> <br> call VecCreate(PETSC_COMM_WORLD,b,ierr)<br> call VecSetSizes(b, PETSC_DECIDE, M,ierr)<br> call VecSetFromOptions(b,ierr)<br> <br> do I=1,M<br> v = 0.5<br> call VecSetValue(b,I-1,v, INSERT_VALUES,ierr)<br> end do<br> <br> call VecAssemblyBegin(b,ierr)<br> call VecAssemblyEnd(b,ierr)<br> </font></div><div><font face="monospace"> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! FIRST KSP/PC SETUP<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> <br> WRITE(*,*) "KSP/PC first setup"<br> <br> call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)<br> call KSPSetOperators(ksp, A, A, ierr)<br> call KSPSetType(ksp, 'preonly', ierr)<br> call KSPGetPC(ksp, pc, ierr)<br> call KSPSetUp(ksp, ierr)<br> call PCSetType(pc, PCGASM, ierr)<br> </font></div><div><font face="monospace"> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! GASM, SETTING SUBDOMAINS<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> <br> WRITE(*,*) "Setting GASM subdomains"<br> <br> ! Let's create the subdomain IS and inflated_IS<br> ! They are equal if no overlap is present<br> ! They are 1: 0,1,8,9<br> ! 2: 10,...,15,18,...,23,...,58,...,63 <br> <br> indices_first_domain = [0,1,8,9] ! corresponds to diag(1,1)<br> do I=0,5<br> do J=0,5<br> indices_second_domain(I*6+1+J) = 18 + J + 8*I ! corresponds to diag(2,2,...,2)<br> !WRITE(*,*) I*6+1+J, 18 + J + 8*I<br> end do<br> end do<br> <br> ! Convert into IS<br> ISLen = 4<br> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain,<br> & PETSC_COPY_VALUES, subdomains_IS(1), ierr)<br> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain,<br> & PETSC_COPY_VALUES, inflated_IS(1), ierr)<br> ISLen = 36<br> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain,<br> & PETSC_COPY_VALUES, subdomains_IS(2), ierr)<br> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain,<br> & PETSC_COPY_VALUES, inflated_IS(2), ierr)<br><br> NSub = 2<br> call PCGASMSetSubdomains(pc,NSub,<br> & subdomains_IS,inflated_IS,ierr)<br> call PCGASMDestroySubdomains(NSub,<br> & subdomains_IS,inflated_IS,ierr)<br> <br> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! GASM: SET SUBSOLVERS<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! <br> <br> WRITE(*,*) "Setting subsolvers for GASM"<br> <br> call PCSetUp(pc, ierr) ! should I add this?<br> <br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,<br> & "-sub_pc_type", "lu", ierr)<br> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,<br> & "-sub_ksp_type", "preonly", ierr)<br> <br> call KSPSetFromOptions(ksp, ierr)<br> call PCSetFromOptions(pc, ierr)<br> </font></div><div><font face="monospace"> <br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br> ! DUMMY SOLUTION: DID IT WORK?<br> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! <br> <br> WRITE(*,*) "Solve"<br> <br> call VecDuplicate(b,x,ierr)<br> call KSPSolve(ksp,b,x,ierr)<br> <br> call MatDestroy(A, ierr)<br> call KSPDestroy(ksp, ierr)<br> call PetscFinalize(ierr)</font><br></div><div><br></div><div>This code is failing in multiple points. At <font face="monospace">call PCSetUp(pc, ierr)</font> it produces:</div><div><br></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><i><font color="#ff0000">[0]PETSC ERROR: Argument out of range</font></i></div><div><i><font color="#ff0000">[0]PETSC ERROR: Scatter indices in ix are out of range</font></i></div><div><i><font color="#ff0000">...</font></i></div><div><i><font color="#ff0000">[0]PETSC ERROR: #1 VecScatterCreate() at ***\src\vec\is\sf\INTERF~1\vscat.c:736</font></i></div><div><i><font color="#ff0000">[0]PETSC ERROR: #2 PCSetUp_GASM() at ***\src\ksp\pc\impls\gasm\gasm.c:433</font></i></div><div><i><font color="#ff0000">[0]PETSC ERROR: #3 PCSetUp() at ***\src\ksp\pc\INTERF~1\precon.c:994</font></i></div><div><i><br></i></div></blockquote><div>And at <span style="font-family:monospace">call KSPSolve(ksp,b,x,ierr) </span><font face="arial, sans-serif">it produces:</font></div><div><font face="arial, sans-serif"><br></font></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><i><font color="#ff0000">forrtl: severe (157): Program Exception - access violation</font></i></div></blockquote><div><br></div><div>The index sets are setup coherently with the outputs of e.g. <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/output/ex71f_1.out" target="_blank">https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/output/ex71f_1.out</a>: in particular each element of the matrix A corresponds to a number from 0 to 63.</div></div></blockquote><div><br></div><div>This is not correct, I believe. The indices are row/col indices, not indices into dense blocks, so for</div><div>your example, they are all in [0, 8].</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Note that each submatrix does not represent some physical subdomain, the subdivision is just at the algebraic level.</div><div>I thus have the following questions:</div><div><ul><li>is this the correct way of creating the IS objects, given my objective at the beginning of the email? Is the ordering correct?</li><li>what am I doing wrong that is generating the above errors?</li></ul><div>Thanks for the patience and the time.</div><div>Best, <br>Leonardo</div></div><div class="gmail_quote"><div dir="ltr" class="gmail_attr"><br></div><div dir="ltr" class="gmail_attr">Il giorno ven 5 mag 2023 alle ore 18:43 Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div><br></div> Added in <b style="color:rgb(200,20,201);font-family:Menlo;font-size:14px">barry/2023-05-04/add-pcgasm-set-subdomains </b>see also <a href="https://gitlab.com/petsc/petsc/-/merge_requests/6419" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/6419</a><div><br></div><div> Barry</div><div><br><div><br><blockquote type="cite"><div>On May 4, 2023, at 11:23 AM, LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</div><br><div><div dir="ltr">Thank you for the help. <div>Adding to my example:<div><i> call PCGASMSetSubdomains(pc,NSub, subdomains_IS, inflated_IS,ierr)<br> call PCGASMDestroySubdomains(NSub,subdomains_IS,inflated_IS,ierr)<br></i></div><div>results in:<br></div><div><div><i> Error LNK2019 unresolved external symbol PCGASMDESTROYSUBDOMAINS referenced in function ... <br></i></div><div><i> Error LNK2019 unresolved external symbol PCGASMSETSUBDOMAINS referenced in function ... <br></i></div><div>I'm not sure if the interfaces are missing or if I have a compilation problem.</div></div><div>Thank you again. </div><div>Best,<br>Leonardo</div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div><br></div> Thank you for the test code. I have a fix in the branch <span style="color:rgb(51,50,56);font-family:"GitLab Sans",-apple-system,BlinkMacSystemFont,"Segoe UI",Roboto,"Noto Sans",Ubuntu,Cantarell,"Helvetica Neue",sans-serif,"Apple Color Emoji","Segoe UI Emoji","Segoe UI Symbol","Noto Color Emoji";font-size:14px;background-color:rgb(251,250,253)"> </span><a href="https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d" style="box-sizing:border-box;font-variant-ligatures:none;color:rgb(31,117,203);text-decoration:none;font-size:13.3px" target="_blank">barry/2023-04-29/fix-pcasmcreatesubdomains2d</a> with merge request <a href="https://gitlab.com/petsc/petsc/-/merge_requests/6394" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/6394</a><div><br></div><div> The functions did not have proper Fortran stubs and interfaces so I had to provide them manually in the new branch.<br><div><br></div><div> Use</div><div><br></div><div> git fetch</div><div> git checkout <a href="https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d" style="box-sizing:border-box;font-variant-ligatures:none;color:rgb(31,117,203);text-decoration:none;font-size:13.3px" target="_blank">barry/2023-04-29/fix-pcasmcreatesubdomains2d</a></div><div> ./configure etc</div><div><br></div><div> Your now working test code is in src/ksp/ksp/tests/ex71f.F90 I had to change things slightly and I updated the error handling for the latest version.</div><div><br></div><div> Please let us know if you have any later questions.</div><div><br></div><div> Barry</div><div><br></div><div><br></div><div><br></div><div><div><br><blockquote type="cite"><div>On Apr 28, 2023, at 12:07 PM, LEONARDO MUTTI <<a href="mailto:leonardo.mutti01@universitadipavia.it" target="_blank">leonardo.mutti01@universitadipavia.it</a>> wrote:</div><br><div><div dir="ltr">Hello. I am having a hard time understanding the index sets to feed PCGASMSetSubdomains, and I am working in Fortran (as a PETSc novice). To get more intuition on how the IS objects behave I tried the following minimal (non) working example, which should tile a 16x16 matrix into 16 square, non-overlapping submatrices:<div><br></div><div>#include <petsc/finclude/petscmat.h><br></div><div>#include <petsc/finclude/petscksp.h><br>#include <petsc/finclude/petscpc.h><br> USE petscmat<br> USE petscksp<br> USE petscpc<br> <br> Mat :: A<br> PetscInt :: M, NSubx, dof, overlap, NSub<br> INTEGER :: I,J<br> PetscErrorCode :: ierr<br> PetscScalar :: v<br> KSP :: ksp<br> PC :: pc<br> IS :: subdomains_IS, inflated_IS<br> <br> call PetscInitialize(PETSC_NULL_CHARACTER , ierr)</div><div> <br>!-----Create a dummy matrix<br> M = 16<br> call MatCreateAIJ(MPI_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,<br> & M, M,<br> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,<br> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,<br> & A, ierr)<br> <br> DO I=1,M<br> DO J=1,M<br> v = I*J<br> CALL MatSetValue (A,I-1,J-1,v,<br> & INSERT_VALUES , ierr)<br> END DO<br> END DO<br> <br> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY , ierr)<br> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY , ierr)<br> </div><div>!-----Create KSP and PC<br> call KSPCreate(PETSC_COMM_WORLD,ksp, ierr)<br> call KSPSetOperators(ksp,A,A, ierr)<br> call KSPSetType(ksp,"bcgs",ierr)<br> call KSPGetPC(ksp,pc,ierr)<br> call KSPSetUp(ksp, ierr)<br> call PCSetType(pc,PCGASM, ierr)<br> call PCSetUp(pc , ierr)<br> <br>!-----GASM setup<br> NSubx = 4<br> dof = 1<br> overlap = 0<br> <br> call PCGASMCreateSubdomains2D(pc, <br> & M, M,<br> & NSubx, NSubx,<br> & dof, overlap,<br> & NSub, subdomains_IS, inflated_IS, ierr)<br><br> call ISView(subdomains_IS, PETSC_VIEWER_STDOUT_WORLD, ierr)<br> <br> call KSPDestroy(ksp, ierr) <br> call PetscFinalize(ierr)</div><div><br></div><div>Running this on one processor, I get NSub = 4.</div><div>If PCASM and PCASMCreateSubdomains2D are used instead, I get NSub = 16 as expected.</div><div>Moreover, I get in the end "forrtl: severe (157): Program Exception - access violation". So:</div><div>1) why do I get two different results with ASM, and GASM?</div><div>2) why do I get access violation and how can I solve this?</div><div>In fact, in C, subdomains_IS, inflated_IS should pointers to IS objects. As I see on the Fortran interface, the arguments to PCGASMCreateSubdomains2D are IS objects:<br></div><div><br></div><div> subroutine PCGASMCreateSubdomains2D(a,b,c,d,e,f,g,h,i,j,z)<br> import tPC,tIS<br> PC a ! PC<br> PetscInt b ! PetscInt<br> PetscInt c ! PetscInt<br> PetscInt d ! PetscInt<br> PetscInt e ! PetscInt<br> PetscInt f ! PetscInt<br> PetscInt g ! PetscInt<br> PetscInt h ! PetscInt<br> IS i ! IS<br> IS j ! IS<br> PetscErrorCode z<br> end subroutine PCGASMCreateSubdomains2D<br></div><div>Thus:</div><div>3) what should be inside e.g., subdomains_IS? I expect it to contain, for every created subdomain, the list of rows and columns defining the subblock in the matrix, am I right?</div><div><br></div><div>Context: I have a block-tridiagonal system arising from space-time finite elements, and I want to solve it with GMRES+PCGASM preconditioner, where each overlapping submatrix is on the diagonal and of size 3x3 blocks (and spanning multiple processes). This is PETSc 3.17.1 on Windows.</div><div><br></div><div>Thanks in advance,</div><div>Leonardo</div></div>
</div></blockquote></div><br></div></div></div></blockquote></div>
</div></blockquote></div><br></div></div></blockquote></div></div>
</blockquote></div><br clear="all"><div><br></div><span>-- </span><br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote></div>
</div></div>
</blockquote></div>
</div></blockquote></div><br></div></div></body></html>